Нестационарный массообмен в вакуумном ротационно-пленочном аппарате при влагоудалении из фосфолипидных эмульсий Текст научной статьи по специальности «Математика»

ВестникВГУИТ, № 1, 2012

УДК 519.8: 665.37.047.79 Докторант C. Алтайулы, профессор С. Т. Антипов,

(Воронеж. гос. ун-т. инж. технол) кафедра машин и аппаратов пищевых производств, тел. (473) 255-38-96

доцент И.О. Павлов

(Воронеж. гос. ун-т инж. технол.) кафедра информационных технологий моделирования и управления, тел. (473) 255-25-50

Нестационарный массообмен в вакуумном ротационно-пленочном аппарате при влагоудалении из фосфолипидных эмульсий

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

The question of creating a computer model designed for the investigation of unsteady mass transfer processes in the dripping of the phospholipid emulsion of sunflower oil in a vacuum cylindrical rotary-film apparatus. The implementation of the mathematical model performance using the finite element method in the Maple.

Ключевые слова: математическое моделирование, метод конечных элементов, тепломассообмен, фосфолипидные эмульсии, ротационно-пленочный аппарат.

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

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

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

© Алтайулы С., Антипов С.Т., Павлов И.О., 2012

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

Нагреваемый участок аппарата можно представить в виде тонкослойной кольцевой фигуры радиусом Я, длиной Ь и толщиной 8 = Я - Яс. Температура внутренней стенки корпуса аппарата на этом участке нагрева поддерживается постоянной и на поверхности пленки происходит испарение под вакуумом.

Выбираем цилиндрическую систему координат, при этом ось г совпадает с осью орошаемого тонкопленочного вращающегося кольцевого слоя в ротационно-пленочном аппарате (г = 0 соответствует положению обрабатываемого продукта в начале рабочего нагреваемого участка).

Уравнение массопроводности для нестационарного процесса массопереноса в цилиндрической системе координат (г', г') имеет следующий вид:

dc = D,

dt

я2 Л д c

1 д ( , дс Л r' дr' і дr' І дz'2

+ W (t), (1)

г ' е[0, Ь], г' е[Яc, Я], t е [0, tk ], где Ят - коэффициент массопроводности; Ж ^) - внутренний источник (сток) распределяемого вещества в единицу времени; йс / й -полная производная концентрации по времени

йс дс дс дс ...

— =—+ уг—т + V г—т; (2)

й д t д г д г

где vr, vz - компоненты вектора скорости в направлении осей r' и z', которые могут быть определены из уравнений Навье - Стокса; tk - время протекания процесса.

Для решения уравнения (2) необходимо задать начальное условие

t = 0,с = Со(z', r',0) , а также граничные условия первого, второго и третьего рода.

Граничные условия первого рода. Задается распределение концентрации вещества у поверхности тела сп (t)

c(z' , r' , t) = сп (t),

например, на поверхности ^1 (рис. 1) задана постоянная концентрация с0 = const:

c(z' , r' , t) = с0.

Граничные условия второго рода.

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

- DmP

С дс Л

дп

= 0

(3)

дс дс , д ,

где -----= -----1г +---- производная от кон-

дп дГ ’ дг'

центрации по нормали п около поверхности тела; 1Г, 1г - направляющие косинусы единичного нормального вектора к поверхности тела; р - плотность распределяемого вещества. Условие (3) используется для задания условия непроницаемости на границе 54 :

дс

= 0.

(4)

S 4 =0

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

+ в(п - 0 = °;

(5)

закон (например, закон Нернста, Щукарева и др.) [4]. Граничное условие третьего рода, заданное на поверхности «2, можно записать в следующем виде:

- я д

т л

. дп ,

\ /и

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

Введем новые переменные

С = с/с0, г = г'/Ьх, г = г'/,т = t/ 1к (6) и учитывая, что внутренний источник отсутствует Ж^) = 0 , получим уравнение (1) в безразмерном виде

йС Г 1 д { дС\ д2С 1

d т

=K

-I r

I + -

r дr ^ дr ) дz

2

z є [0, zk], r є[гс, 1], тє [0,1], с начальным условием C (z, r,0) = 1 и граничными условиями

C ( z, r,T)l z=0 = 1

д^ I = Ki •

дr Ir=l m ’

_dC_\

дr

(7)

(8)

(9)

(10) (11)

= Б1т (С|г=гс - О, ), где К1 = 1/ (Яе- 8с) - комплексный критерий; Re=VxL/v - число Рейнольдса; 8с = у/ Ят -число Шмидта; Б1т = в / Ят - критерий Био (массообменный); К1т = (/) Ьгл /фт рс0)

- критерий Кирпичева (массообменный); 2(г) = Ж^) tk / с0 - источник (сток) массы в безразмерном виде; V - коэффициент кинематической вязкости; Сж = сх / с0 . Полная производная концентрации по времени определяется соотношением

dC

дC дC дC

' + Ur--------+ Uz

И д ' д - д ’ (12)

ат дт дг дг

где иг = уг / , иг = - масштаб-

ная длина области и базовая скорость системы. За масштабную длину удобно принять наружный радиус пленки = Я, тогда гс = Яс / Я,

гк =1 / Я^ =1 / Ь .

Дифференциальное уравнение массопроводности (7)-(12) будем решать методом конечных элементов (МКЭ) [5, 6, 7, 8, 9], который основан на идее аппроксимации непрерывной функции (в данном случае - концентрации) дискретной моделью, строящейся на

ВестнщВГУИТ, № 1, 2012

множестве кусочно-непрерывных функций, определенных на конечном числе подобластей, называемых элементами [7].

Решение дифференциального уравнения теплопроводности (3) с граничными условиями (4)...(9) заменим другой задачей - найти распределение концентрации С (г, г), которое минимизирует следующий функционал:

з -12 2 Vа

К

дг

+

Віт

дг

+2сГ аС

I Ит

г dV +

+ К1т | гС + —т |г (С - Сх)2 , (13)

«2 2 «3

где V - объем исследуемой области; « - площадь поверхности.

Введем матрицы

[ё]Т = [дС / Иг дС / Иг] и [р] =

гК2 0

0 гК,

.(14)

С учетом (14) функционал (13) принимает вид

З = 1ё ]Т [Р][ё] dV + V С^г dV +

+ Кіт 1 гСИБ + В2т 1 г (С - О2 ИЗ . (15)

^ 2 Зз

Вся область разбивается на ЫЕ число конечных элементов, тогда функционал (15) можно представить в следующем виде:

ЫЕ ( 2

з -Е Зе 2, (16)

е-1

где З(е) - функционал, записанный для конечного е -элемента, выражение которого имеет такой же вид, как и для функционала (15):

з - Ті 2 [ <е> Г [о « ] <е> ] -1 с <е> ^ +

е-1 Vі

+ Кіт 1 С^ИЗ + ^ 1 г(с(е) -С2} ИЗ , (17)

С(е) -

концентрация в е -м конечном эле-

- 0. (18)

где менте.

Условие минимума функционала (17) требует выполнения соотношения

дЗ д хе (е) ш д/?)

д{0} д{0} е=1 е=1д{С}

Частные производные дЗ(е) / д{~} в (18) не могут быть вычислены пока интегралы в (17) не будут выражены через узловые значения концентрации распределяемого вещества {~}.

Учитывая соотношение интерполяционного полинома, записанного через функции формы

ГС1

С(е) - [(е)] {}- [,[2е)...Ыре)

С2

Сг

, (19)

В

(е)

где р - число узлов элемента, можно вычислить величину [ g(е)] (17):

[ g(е)] = |в(е)|{~}, (20)

где | В(е) | содержит информацию, связанную с производными от функций формы

дЫ<е)/ дг д^2е)/ дг ••• дЫ(р)/ дг

ды1е) /дг ды2е)/дг ••• дЫ^ /дг

Использование формул (19) и (20) позволяет записать интегралы по элементам в выражении функционала (17):

з (е) = | -2 {~} [в(е)]г [я(е)][в(е)]{~} (IV+

V 2

!Ме), д С(е) д С(е) д С (е). _

+ | С (——------+ Ыг —“----+ и г —“--) г (V +

V

+

дт д г

1 Кіт [ы(е)]{~}гИБ +

д г

З2

+

1 В2т {~Т [ы (е) Т[ы (е) ]{~~}гИЗ-

-1 В^С„[ы(е) ]с}гИ8 +1 Віт

Зз

2

С2 г ИЗ. (21)

Из условия минимального значения для функционала (21) получим систему линейных дифференциальных уравнений вида

т(е) ||Т} + [к(е) 1~}={(е)}, (22)

где [т(е)]= | г {[Ы ] [Ы ]1 (V ;

V(е) к ;

1к (е)]=|к1е)]+|к2е)]+к3е)]-

|к('*]= ( [В('>]г[D<^)][B^‘>]c)dV;

V(е)

к2е)]= |г[Ы] [Ы] Б1т (А;

[Ы ] иг + и_ ^1- г (V :

I-

V

( " дЫ" " дЫ" \

и,„ + и

г дг г дг

{р(е 2}-{р1(є2}-{р2є 2}+{р3е 2};

з

З

4

2

{p(e)}= 1 r[N]TQdV ;

v(e)

При этом производные функций форм по координатам определяются как

{p2e)(T)}= I r [N]T Kim dA .

s(e)

Общая система уравнений формируется из систем линейных дифференциальных уравнений конечных элементов (22):

[M { }+[^ ]{c}={p}. (23)

Систему (6) решаем, используя схему Кранка - Николсона [9]. Предполагаем, что

{г+дЛ-{сЛ = ^({+Дг}-{}), (24)

где Ат - шаг интегрирования. Из (24) выражаем скорость изменения концентрации в момент времени т + Дт в виде следующего уравнения:

{+ДГ}={}+Д-({СГ+ДГНСЛ). (25)

Подставив выражение (25) в (23), получаем конечно-разностное уравнение для определения искомой зависимости

—м+K Дт

{Ст+Дт}=[М ](ДтСтМСт})+К },(2б)

где {/ср} - вектор-столбец правой части уравнения (26) в момент времени Т + А т / 2 . Алгоритм этой схемы состоит в последовательном решении уравнения (26).

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

Для аппроксимации концентрации С ( г) во всей области конечного элемента используем функции форм следующего вида:

с (p, 5) = У Ni (p, і)сі = [n (p, і){^ ^

где Nt (p, 5) - функции форм;

Ni ( 5 ) = 4 (l- p)(l - 5);

N 2 ( 5) = 4(l+p)(l - 5); N3 (5 )=1 (l+p X1+5); N4 ( 5) = 4 (l - P № + 5).

dN dN

дт dN = [J і-1 д p dN , ( = (...,4

д z д s

где [j ] - матрица Якоби следующего вида: [j ] =

4 N 4 dNі

y^L r У~Т*

і=i dp і=i dp

4 dN, 4 dNі

У і r у d 4 Удz і =1 д

где г, - глобальные координаты узлов ко-

нечных элементов. Элементарный объем СУ вычисляется как

СУ = ёе1:[/] ёр Ся , где А(P, я) - толщина элемента;

А(р, Я)=[Ж(р, 5)]АА и =[а! а2 А А4] .

Элементарная площадь поверхности ёА вычисляется различно для сторон КЭ: для сторон 1 и 2

СА = ^(дг/ др)2 + (д г/ др)2ёр ;

для сторон 3 и 4

СА ^ ^(дг7д^)2+"(дй7д^)2"dq .

Матрицу вида [т(е)] назовем "матрицей массоемкости” (по аналогии с названием этой матрицы “матрица теплоемкости” при расчете процессов теплообмена [5]). Общая система уравнений формируется из систем линейных дифференциальных уравнений.

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

Рис. 1. Изменение концентрации распределяемого вещества по толщине пленки и длине аппарата: 2, г - безразмерные пространственные координаты; wcp - средняя скорость движения слоя, м/с

С, кг/кг 0.5

0.4

0.3

0.2

0.1

0 0.5 1 1.5 2 2.5 3 3.5

------► г

’Л'ср

Рис. 2. Изменение влажности фосфолипидной эмульсии на поверхности массообмена по длине аппарата: 2 - безразмерная пространственная координата; wcp - средняя скорость движения слоя, м/с

Использование математического пакета символьной математики Мар1е позволяет эффективно реализовать алгоритмы метода конечных элементов для решения задач тепло- и массопереноса. Конечно-элементный анализ нестационарного массообмена в вакуумном ротационно-пленочном аппарате при влаго-удалении из фосфолипидных эмульсий позволяет выполнять проектные процедуры

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

ЛИТЕРАТУРА

1. Пат. № 99987 РФ на полезную модель МПК ВО1Б 1/22. Цилиндрический ротационно-пленочный аппарат [Текст] / Алтайулы С., Антипов С.Т., Шахов С.В.; заявитель и патентообладатель Воронеж. гос. технол. акад. -№ 2010110753/05; заявл. 22.02.2010.; Опубл. 10.12.2010, Бюл. №34 - 4 с.

2. Кафаров, В. В. Основы массопередачи [Текст] / В. В. Кафаров. - М: Высш. школа, 1979. - 439 с.

3. Лыков, А. В.Теория тепло- и массопереноса [Текст] / А. В. Лыков, Ю. А. Михайлов.

- М.: Госэнегоиздат, 1963. - 563 с.

4. Рудобашта, С. П. Массоперенос в системах с твердой фазой [Текст] / С. П. Рудобашта.

- М.: Химия, 1980. - 248 с.

5. Аладьев, В. З. Маріє 6. Решение математических, статистических и физико-технических задач [Текст] / В. З. Аладьев, М. А. Богчавичус. -М.: Лаборатория базовых знаний, 2001. - 824 с.

6. Зенкевич, О. Метод конечных элементов в технике [Текст] / О. Зенкевич. - М.: Мир, 1975. - 541 с.

7. Сигерлинд, Л. Применение метода конечных элементов [Текст] / Л. Сигерлинд. -М.: Мир,1979. - 392 с.

8. Скульский, О.И. Механика аномально вязких жидкостей [Текст] / О. И. Скульский,

С. Н. Аристов. - Москва-Ижевск: НИЦ “Регулярная и хаотическая динамика”, 2003. - 156 с.

9. Митчелл, Э. Метод конечных элементов для уравнений с частными производными [Текст] / Э. Митчелл, Р. Уэйт. - М.: Мир, 1981.

- 216 с.