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

УДК 664.002.05

Математическая модель массопереноса при экстрагировании слоя

Кошевой Е.П., Меретуков З.А., Косачев В.С., Михневич А.Н.

[email protected]

Майкопский государственный технологический университет, Кубанский государственный технологический университет

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

Ключевые слова: растительный материал, массообмен, экстрагирование, математическое моделирование.

Mathematical model of mass transfer at extraction process of a

layer

Koshevoy E.P., Meretukov Z.A., Kosachev V.S., Mihnevich A.N.

[email protected]

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 offinal differences.

Keywords: vegetative material, mass transfer, extraction, mathematical modelling.

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

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

Главные предположения о модели следующие:

- система является изотермической и изобарической;

- на входе в слой физические свойства растворителя постоянные;

- радиальные градиенты концентрации в жидкой фазе отсутствуют;

- перемешивание жидкой фазы имеет место только в осевом направлении;

- экстракт принят как единственный компонент;

- концентрация экстрактивных веществ в твердых частицах изменяется только в радиальном направлении и не зависит угла направления радиуса.

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

Уравнение переноса в твердой фазе определяется уравнением вида

дт

2 L

pep R

£

дС,

(1)

-5

где Cs - концентрация экстрактивных веществ в твердой фазе, кг/м ; т безразмерное время (т = lJ0t/Ls); U0 - скорость жидкости в расчете на незаполненное сечение экстрактора, м/с; t - время, с; L - длина слоя, м; в порозность слоя; Рер - число Пекле частицы твердой фазы (Рер = U0dp/Dms)', dp - диаметр частицы, м; Dm - коэффициент внутренней диффузии в твердой фазе, м/с; R - радиус частицы, м; с, - безразмерный радиус частицы (t = r/R). Преобразуем уравнение (1) в виде суммы производных

8CS _ 2 L d2Cs 4 L dCs дт ~ Pep R 84' - -p

de

(2) При граничных условиях г>0, £=1, -

(3)

öcf

Pep R

= Bi • С, ~Cf .

,cfi=kp.c:

Уравнение массопереноса в жидкой фазе определяется зависимостью вида

dCj_

дт

1 d2cf de

Реь gZ'

dZ

Bi

e-R Per

'ft .

C/"

Peb dZ

= 0

(4)

при т> 0, 2=0,

(5)

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

С -С

]+1 }

(7)

Э^2

С -2-С +С

7-1,У+1 ^ /+1,у+1

д#2

Для аппроксимации первой производной по координате будем использовать центральную аппроксимацию для внутренних точек сетки в твердой фазе

2-А^

(8)

Для аппроксимации первой производной на границах сетки будем использовать левую или правую аппроксимации соответственно

(9)

С -С

С1,]+1 С1+1,;+1

(10)

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

Обычно безразмерная координата аппроксимируется выражением

£ = 1-А^

(П)

Однако в нашем случае необходимо произвести подбор разбиения координатной оси исходя из следующих условий

1. Шаг А£, должен обеспечивать устойчивость алгоритма решения линейной системы уравнений.

2. Граница сопряжения твердой и жидкой фазы должна совпадать с узлом сеточной схемы.

Для типичных параметров экстракции этим условиям соответствует значение Л£,=2/10. Узлы разбиения берутся по формуле (11) для 1=0...5. В этом случае узлы сетки представляют собой зависимость табличного вида

Таблица 1 Узлы сетки при сеточной аппроксимации

1 -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 Ъ

0 1 2 3 4 5 6

где {- номер узла сетки, с, - значение текущей безразмерной координаты в твердой фазе, Ъ - значение текущей безразмерной координаты в жидкой фазе.

Точки для 1=(-1) и 1=12 используются для аппроксимации условий симметрии на границах соответствующих фаз.

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

2 Ь С,-, ,-+, - 2-С. . , +С +] . , 4 Ь 1

С -С

С -С

Сг-1, ]+1 Сг+1, ]+1

Рер Я

Л#2

Рер Я

2 ■

(12)

Данная схема содержит член который делает веса этой схемы нерегулярными. Следовательно, каждый из этих членов должен рассчитываться индивидуально по каждому узлу.

Начнем построение сетки и вычисление её весов с нулевого узла. Характерными особенностями этого узла являются сингулярность и условие симметрии.

Учитывая, что функция концентрации в этой точке ограничена, то значение выражения в этой точке имеет предел равный 0, т.е. имеет место, следующее выражение

Шп

= 0

(13)

Следовательно, при сеточной аппроксимации для 1=0 член, содержащий первую производную по концентрации, отсутствует, тогда

2 1Си+1-2-С,

С,

О.у+1

■с,

0,1

0,у+1

■С

1.У+1

Рер Я

(14)

Условие симметрии в этом узле заключается в равенстве концентраций слева и справа от этого узла. Поэтому окончательно имеем

2 Ь 2- С1;+1 — 2 •

С -С

0,7+1 0,7

РерЯ

А^

(15)

Для реализации неявной схемы необходимо произвести разделение временных слоёв относительно знака равенства 0 и ]+1 члены) и преобразовать данное уравнение в схему с весами при этих членах. Проводя эквивалентные алгебраические преобразования, имеем:

Со,7 =

Рер -Я- АВ,2 + 4-Ь-Ат Рер -Я-А^2

' Со,7+1 +

■4-Ь-А т

Л

Я- АВ,2

■сх

У+1

о

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

Для промежуточных узлов в пределах от 1= 1,___,4 алгебраические

уравнения неявного вида для этих узлов получаются подстановкой соответствующего индекса (номера узла по координате) в опорное уравнение (12).

Для первого узла (¿=1) имеем

С -С

^1,7+1

2 Ь Со,^ • + С2>;-+1

РврЯ

ц2

+ ■

4 Ь 1

Рер Я г- Ц

2 ^0,у+1

■С,

У+1

(17)

и после аналогичных преобразований получаем следующее трехчленное разложение

си =

( -6-Ь-Ат Л

■ С0>;+1 +

(Ре -Р-Д^2 +8-1-ДгЛ

Ре -Р-Д£2

'Си+1 +

Г 2-Ь-Ат Л

Ре -Р-Д£2

V р

■С

2,7+1

(18)

В остальных внутренних узлах сетки по твердой фазе будем иметь аналогичные уравнения следующего вида. Для второго узла (1=2) имеем

-Ъ-Ь-Ат

Рер - Р-Д^2

•Си+1 +

ч р (19)

Для третьего узла (1=3) имеем

Рер -Р-Д^2 +4-Р-Дт Рер - Р-Д^2

•С2,7+1 +

Ь-Ат

Рер - Р-Д^2

ч р

С

3,7+1

^ =

Г -%-1-Ат А

3-Рер-Я-А%2

Г Рер-Я-А^2 +4-1-АтЛ

Рер-Я-А%2

А-Ь-Ат

З-Ре -Я-Ад

С

4,7+1

(20)

Для четвертого узла (1=4) имеем

( -5-Ь-Ат Л 2■Рер-Я-А^2

Рер -Я-А^2 +4-Ь-Ат Рер-Я-А%2

' ^4,7+1 +

З-Ь-Ат

2-Ре -Я-Ад

■С

5,7+1

(21)

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

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

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

соприкасающейся с жидкой фазой. В дифференциальном виде на границе твердой фазы должны выполняться следующие соотношения

дС.

а?

= Ш ■ {'/V "С/.

С/в ~ кр '

При г > О, £=1,

(22) и (23)

где В1 - число Био (В1 = к^Е/От); к] - коэффициент массопередачи, м/с; кр -объемный коэффициент распределения; С^з - концентрация экстрактивных

-5

веществ в жидкой фазе на поверхности частицы, кг/м ; - концентрация экстрактивных веществ в жидкой фазе, кг/м3; С+ - равновесная концентрация

-5

экстрактивных веществ в твердой фазе на поверхности частицы, кг/м .

Учитывая необходимость сохранения диагональной структуры весов матрицы, объединим эти уравнения следующей разностной схемой

1С

-Г

4,у+1 6,7+1

2

Р 5,7+1 6,7+1 .

(24)

В опорном уравнении (12) присутствует член, содержащий значение первой производной, который для данного граничного условия (24) принимает следующий вид

-4,7+1 Сб,7+1

'=-2-А--С -С:

р 5,7+1 6,7+1

с.

(25)

В результате на границе раздела фаз со стороны твердой фазы уравнение (12) принимает следующий вид

С -С

2 Ь Сг_1 ,+1 - 2 • С, ,+1 + С.

7+1 7+1,7+1

РерЯ

Л<?2

+ ■

4 £ 1

Ре Л г•

^ -С

Р 5,7+1 6,7+1 ^

(26)

Подставив номер граничного узла (1=5) в это уравнение получаем С-5.7+1 — > 2 Ь С-4. ,+1 — 2 • С,+ С6;+] > 8 Ь ■ Вт * _ _

■• -с,-а

Аг

РерЯ

5-Рер Я-А^

*р 5,7+1 6,7+1 .

(27)

Или в виде удобном для использования в решении

-2-Ь-Ат Ре-Я- Д£2

5 ■ Рер ■ Я ■ А^2 + 20- Ь- Ат + Ь- Вг - кр ■ Ат ■ Рер -Я-А?

5 • Рер • Я • А^2

(28)

Следующие узлы разностной схемы принадлежат жидкой фазе (1=6,...,11), для которой на границе с твердой фазой (1=6) выполняется дифференциальное уравнение (4)

1 дС

при т>0, С,—---^ = 0

^ 1 Реь д2

(29)

где Ъ - безразмерная осевая координата по слою, 7/Ь; z расстояние измеренное от входного отверстия слоя, м; Реь - число Пекле жидкой фазы в слое частиц (Реь = и^Ю^е)', /Л. коэффициент осевой дисперсии жидкой фазы в слое частиц, м2/с.

Преобразуя значение производной в уравнении (29) в разностную схему для этого узла получаем уравнение граничного условия на разделе фаз со стороны жидкой фазы

1 С5,;+1- с 1,]+\ ^

с -

б';+1 Реъ 2А7

(30)

Уравнение массопереноса в жидкой фазе имеет вид

дСг _ 1 д2сг дСг £ _с -

Вт ~ Реь дI2 32 + е-Я Рер * ^

(31)

Преобразуем данное уравнение (31) в разностное уравнение

Сщ ~ _ 1 ~ ^ ' + ~ ^¿+1,;+! 6 • (~ £ Л, • В! I Г _Г ^

Ат ~Реь к!1 " 2-Аг £-Я-Рер р' 5>}+1~

(32)

При этом на границе значение первой производной входящее в это уравнение может быть заменено подстановкой граничного уравнения (30)

~ о о„ -Р ^,./+1 ■

Дг дг2 и'7Т1 " г-Д-Ре,

(33)

Подставив в него номер граничного узла (1=6) получаем

_ 1 С52• Сб>;.+1 + С76- {-£ф _

Дг дг2 и'7Т1 " г-Д-Ре,

(34)

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

Ь ■ В1 ■ Ре, -е-к ■ Ат ■ Аг2-е-Я-Ре ■ Ат-б-Ь-Вг-Ре, -к ■ Ат ■ А22

С ___р_р_° р_ с .

е-Я-Реъ.Рер. Ь22 >:Л "

Реъе-Я-Ре№2 +2еЯ-Ре Ат + Ре2е ■ Я-Ре АтАг2 + | ■ Ь ■ Вг ■ РеъАт№2 ^ <-

^ о р р о р ^ о ^ С

"' е-Я-Реь-Рер-Аг2 '

+ __С

Реъ-А22 Ъ1+1

Как видно из представленного уравнения (35) сохранена трех диагональная структура, что позволяет получать на каждом временном слое взаимно-согласованное решение по обеим фазам. Внутренние узлы содержат член с концентрацией в граничном узле твердой фазы. Для седьмого узла (¿=7) имеем

~ _ 1 С6и+1 - 2 • С7>;.+1 + - С8>;.+1 ^ 6 ■ (- ■ Ш л _ -

'4,7+1 4,7+1 „

Ре,

А7

2-дг

е-Я-Ре.

(36)

Для восьмого узла (¿=8) имеем

^8,7+1 - О/ 1 _ 2 • С8 + С"9;7+1 Су - С"9;7+1 1 6 • С - £' ; В/

Ре и

(37)

Для девятого узла (1=9) имеем С -С 1 Г -2-Г +Г

9. /+1 9 / -I ^8,7+1 ^

_ 9]

1-А1

С -С

Б-Я-Ре,

I •С -С

5,7+1 8,7+1 .

Ре„

А!

9,7+1 1 10,7+1 8,7+1 10,7+1 6 -{-е^-ВХ л -

2 " 1-А1 + е-Л-Ре. '^'4+1-4+1.

(38)

Для десятого узла (1=10)

0,7+1 ~ ^10,; 1 ^9,7+1 ~ ^ ' 0,7+1 + 1,7+1 ^9,7+1 ~ ^11,7+1 6 • ( - £ • Ь • В!

Рек

М

'^в'С 5./+1 С]

р 5,/+1 10,7+1.

(39)

Проводя эквивалентные алгебраические преобразования этих уравнений к виду удобному для решения, имеем Для седьмого узла ^=7)

С-и =

-б-Ь-Ш-Ре^ - к„ • Ат ■ Д2Г2 + 6-Ь-Ш • Ре, к-£-АтА2

"Ъ р

Ъ р

а-Я-Реь ■Рер ■ А22 Реь - а-Я-Рер ■Ат-А2-2-е-Я-Рер ■ Ат 2-е-Я-Реъ ■Рер ■

-с,

Реъ ■£■ Я- Ре р А22 +2-е-Я-Рер ■ Ат + б-Ь-Вг ■ Реь ■Ат■Аг2 ■ К-е

Реь ■ Я ■ Рер е ^2

-Реъ ■£■ Я-Рер ■Ат-Аг-2-е-Я-Рер Ат •••н 7г

•С

7,7+1

2-Реь ■£■ Я-Рер Д^

(40)

Для восьмого узла ^=8) имеем

С8,7 =

6-Ь-В1-Ре, -к„ Ат-Аг2 +6-Ь-В1-Ре, -к„ -е-Ат-Аг

Ър

Ър

£■ Я-Реь ■Рер -Лг2

Ре, ■£-Я-Реп ■Ат-Аг-2-£-Я-Реп-Ат

-+—-£-2-£--С7/+1+...

2-£-Я-Реь-Ре -А22 ''

2

2

Реь - е- Я- Ре р ■ А22 + 2-£-Я-Рер ■ Ат + 6 ■ Ь ■ Ш ■ Реь ■Ат■Аг2 ■ 1-е

Реь -Я-Рер-£-№2

•С.

8,7+1

-Реь -£-Я-Рер-Ат-К1-2-£-Я-Рер-Ат

2-Реь ■£■ Я-Рер

■С,

9,7+1

(41)

Для девятого узла (1=9) имеем

-6-Ь-Вг-Реь -кр ■ Ат ■ №2 + 6 ■ Ь ■ Ш ■ Ре, -к„ £-АтА2

Ь р

£-Я-Реъ ■Рер ■

Ре, ■Б-Я-Реп -Ат-Аг-2-£-Я-Ре-Ат +_-_I_I__с +

2-£-Я-Реь-Ре -Аг2 ^

Реь - е-Я- Ре р ■ ¿£2 +2-е-Я-Рер ■ Ат + 6 ■ Ь ■ Ш ■ Реь ■Ат■Аг2 ■ К-е

Реь - Я-Рер-е- Дг2

•С,

9,7+1

-Реь 8-Я-Рер ■Ат-Аг-2-£-Я-Рер-Ат

•••Н Гг ^10,7+1

2-Реь ■£■ Я-Рер Д^

(42)

Для десятого узла (1=10).

-б-Ь-Вг -Ре, -к -А т ■ Аг2 +6 ■ Ь ■ Вг ■ Ре, -к е-Ат-А^2 ^ __° р_° р_

1<и ~ £-Я-Реъ-Рер-№2 '

Ре, -е-Я-Ре ■Ат-Аг-2-£-Я-Рев ■ Ат +_-_I_I__С +

2-а-Я-Реь-Ре -Аг2 ^

5,7+1

Реь ■£■ Я- Ре р -Аг2 +2-£-Я-Рер ■ Ат + 6 ■ Ь ■ Вг ■ Реь ■ Ат ■ А22 ■ 1-е

РеЬ-Я-Рер-8-А22

' Сю,7+1 + '''

-Реь-£-Я-Рер-Ат-А1-2-£-Я-Рер-Ат "+ 2-Реь-£-Я-Рер-Аг2 1и+1

(43)

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

С -С

*-1,7+1 ¿+1,7+1

2-А7

2

Используем это условие в уравнении на крайнем узле разностной схемы

Ре„

AZ

2

е-Я-Ре.

р

(45)

которое преобразуем к стандартному виду

-6-Ь ■ В1 ■ Ре, ■ к ■ Ат ■ А22 + 6-Ь-В1 ■ Ре, ■ к е-АтА22 ^ __° р_° р_

1и ~ е-Я-Реъ-Рер-А22 '

2-А т

•С,,

5"7+1 Ре,-А21 1<и+1

Реъ ■ е Я- Ре р -Аг2 +2-е-Я-Рер ■ Ат + б-Ь-Вг ■ Реь ■ Ат ■ А2,2

Реь ■ Я ■ Рер е №2

■■с,.

]+1

(46)

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

с

с,к

с

с

с

С5к

с

с с7,к

с "к

с

с

с

е0

й е1

с2 й 2 е2

С3 йз е3

с4 й 4

й5 е5

сб й 6 е6

Л с7 й е7

/8 с8 й8 е8

/9 с9 й9

/10 /11

"10

й

10

"10 й11

г

г

г

г

Г

г с5,к+1

с

г

Г

г ^9,к+\

Г

г Ми+1

(47)

Как видно из структуры представленной разреженной квадратной матрицы решение данного матричного уравнения не возможно методом прогонки, поэтому использовали метод Гаусса. В этом случае матрица не только не может быть вырожденной, но даже и почти вырожденная матрица не может быть использована. Вырожденность определяется как матрица, у которой детерминант равен нулю. Матрица является почти вырожденной, если имеет большое число обусловленности. Для расчета числа обусловленности с(А) использовали формулу

А.

(48)

где X - собственные значения этой матрицы.

Основным параметром, влияющим на число обусловленности, является шаг сетки по времени. Поэтому исследовали его влияние на это число. В

с

е

4

с

5

е

9

с

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

с 4Q= -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 с.