CC BY
43
УДК 519.688 Дата подачи статьи: 25.02.21
DOI: 10.15827/0236-235X.135.466-476 2021. Т. 34. № 3. С. 466-476
Математическая модель вакуумной сублимационной сушки с неравномерным распределением паров по объему камеры
Е.К. Мохова 1, магистрант, lizi5263@yandex.ru М.Г. Гордиенко 1, д.т.н, доцент, gordienko@muctr.ru
1 Российский химико-технологический университет им. Д.И. Менделеева, г. Москва, 125480, Россия
Вакуумная сублимационная сушка применяется для обезвоживания термолабильных веществ и широко используется в биотехнологии, фармацевтической и пищевой промышленности для получения сухих тканей, препаратов и продуктов без потери их структурной целостности и биологической активности. В настоящее время существует необходимость исследования влияния градиентов давления, возникающих в вакуумной камере, на скорость процесса сушки материалов, находящихся на разном удалении от конденсора. Неравномерный отвод паров растворителя от поверхности материала во время сушки может приводить к возникновению существенных различий в скорости сушки в разных участках камеры.
В данной работе были совместно использованы математическая модель для расчета кинетики процесса вакуумной сублимационной сушки (собственное ПО) и модель на основе вычислительной гидродинамики (Computational Fluid Dynamics, CFD) для прогнозирования давления водяных паров над поверхностью материала (программный пакет Ansys Fluent 17.0). Передача данных между моделями осуществлялась посредством передачи значений массового потока испаряемой влаги, рассчитанного по модели кинетики сушки, через файлы UDF в заданные интервалы времени. Данные UDF-файлов использовались в качестве граничных условий. Полученные распределения парциального давления водяного пара передавались в модель кинетики в качестве граничных условий для следующего цикла итераций. Совместное решение двух математических моделей позволило получить кинетические характеристики процесса сушки для 2 масс. % водного альгината натрия с учетом неравномерного распределения водяных паров по объему камеры. Было показано, что для модели лиофилизатора ScanVac Coolsafe 100-9 при рассмотренных массовых потоках наблюдается быстрое выравнивание профиля давления по высоте рабочей камеры, что обеспечивает однородные условия сушки материала. На основании результатов серии вычислительных экспериментов даны рекомендации по тепловому режиму вакуумной сублимационной сушки материалов медицинского назначения на основе альгинатного гидрогеля.
Ключевые слова: вакуумная сублимационная сушка, тепло- и массоперенос, кинетика сушки, CFD, математическое моделирование.
Разработка математических моделей является одной из актуальных задач в области исследования процесса вакуумной сублимационной сушки. Модели строятся на описании тепло- и массообмена процесса сушки, выраженного с помощью дифференциальных уравнений [1, 2]. В результате решения системы дифференциальных уравнений могут быть получены кинетические характеристики процесса сушки - изменение температуры и влаго-содержания материала. Такие математические модели позволяют, варьируя рабочие условия процесса (температуру нагревания полок и глубину вакуума в камере сублиматора), подобрать режимы ведения процесса, обеспечивающие снижение энергетических затрат и времени сушки [3].
За последние несколько лет появился большой интерес к моделированию промышленного и лабораторного оборудования с использованием вычислительной гидродинамики (Computational Fluid Dynamics, CFD). Это связано с тем, что CFD-моделирование можно использовать для улучшения конструкций, оптимизации геометрии и оценки влияния различных вариантов устройства оборудования на ход процесса [4-6].
Недавние теоретические исследования и моделирование газодинамики показали, что на поверхности сублимации возможны градиенты давления [4, 5]. Давление в камере может существенно повлиять на скорость теплопередачи и, следовательно, на скорость сублимации материала [7]. При достаточно высоком градиенте
давления материалы, находящиеся на разных полках, могут достигать разной степени сухости и в результате повлиять на однородность партии или даже повредить продукт.
Цель данного исследования - разработка и программная реализация математической модели процесса вакуумной сублимационной сушки с учетом неравномерного распределения водяных паров по объему камеры и их влияния на кинетику процесса. Разработанная математическая модель позволяет подбирать технологические режимы ведения процесса для различных материалов, что существенно сокращает время сушки и энергетические затраты.
Объединение математических подходов
Для реализации цели разработанная математическая модель кинетики процесса сублимационной сушки решалась совместно с CFD-моделью распределения водяных паров по объему вакуумной камеры. Данный подход позволил учесть влияние распределения паров рас-
творителя над поверхностью материала на разном удалении от входа в десублиматор на скорость процесса.
Для этого значения изменения температуры и влагосодержания материала, полученные по кинетической модели, записывались в UDF-файлы temperature_flow и mass_flow_rate, которые прикреплялись к соответствующим областям созданной SD-геометрии сублимационной камеры в качестве граничных условий в среде моделирования Ansys Fluent 17.0. На рисунке 1 представлен алгоритм связи двух моделей.
Значения давлений, полученные при CFD-моделировании, фиксировались и затем использовались для расчета кинетических характеристик для нижней и верхней полок отдельно. Такой подход позволил получить модель с учетом неравномерного распределения водяных паров по объему вакуумной камеры, а также прогнозировать давление над полками в рабочей камере и исследовать влияние градиентов давления на скорость процесса сушки материалов, находящихся на разном расстоянии от конденсора.
Рис. 1. Алгоритм связи двух моделей Fig. 1. An algorithm of communication between two models
Математическое моделирование кинетики вакуумной сублимационной сушки
При разработке данной математической модели за основу были приняты статьи [8, 9]. Математическая модель является одномерной, то есть перенос тепла и влаги происходит только вдоль оси Х и строится на уравнениях тепло- и массопереноса. В структуру настоящей математической модели входят описания первого и второго периодов сушки.
Для математического описания первого периода сушки были приняты следующие допущения:
- существует фронт сублимации Хр, который делит образец на замороженную (I) и высушенную (II) области;
- граница раздела (фронт сублимации) в процессе сушки равномерно двигается вниз, пока вся замороженная свободная влага не удалится из материала толщиной Ь;
- на границе раздела фаз выполняется условие: массовый поток через границу в направлении снизу вверх равен массовому потоку через границу в направлении сверху вниз;
- высушенная и замороженная области рассматриваются как гомогенные с точки зрения постоянства теплопроводности, плотности и теплоемкости;
- теплопередача к материалу осуществляется кондуктивно от полки снизу.
Далее представлено математическое описание первого периода сушки.
Для замороженной области (I) уравнение теплопроводности выглядит следующим обра-
. dTi
а
зом: —- = ае1 — +-1—, где Т- температура,
д дх2 ре1Сре1
К; t - время, с; ае1 - температуропроводность, м2/с; X - декартова координата, м; Qv\ - объемная мощность источников теплоты, Вт/м3; ре1 -плотность, кг/м3; Сре1 - удельная теплоемкость, Дж/(кг-К).
Уравнение теплопроводности для высушенной области (II) следующее:
dt
dl
C
= a„,
d( NT )
dXl PbullC,
pell
dx
Ô*
Рш1Сре11
где рьып - насыпная плотность сухого материала, кг/м3; - поток водяного пара в области II, кг/(м2-с).
Приращение координаты фронта сублимации АХп+1 и влагосодержания А У*п+1 в безраз-
мерном виде рассчитывалось для областей (I) и (II) соответственно:
АХ^1 =АЩ'Я+1 = ,'1п+1Бгвл Аг',
АХ ¡"+1 = АЖ*\п+1 = чТ+131еЛ1 Аг', где д* - безразмерный тепловой поток, Вт/м2; Stev - модифицированное число Стефана; Аt -безразмерный шаг по времени; п - позиция во времени.
Расчет смещения фронта сублимации Х*п+1 и изменения влагосодержания Уп+1 проводился также в безразмерном виде соответственно:
X *п+1 = X 'п +АХ *п+1 + АХ ."Г1,
Ж *п+1 = у *п *п+1 + АЖ*п+1.
Перевод влагосодержания в размерный вид осуществлялся по формуле
у = у * (Ур - Ж,) + Ж,,
где V - безразмерное влагосодержание; 1р -начальное влагосодержание материала, кг/кг; Уед - равновесная влажность в первый период сушки, кг/кг.
Начальное условие для замороженной и высушенной областей материала: Т = Тп, 0 < х < Ь, t = 0.
Представим граничные условия для нижней (I) и верхней (II) частей материала: То = Ть х = 0, t > 0, Тп = Ть, х = Ь, t > 0, х = X(t), t > 0, где Т0 - температура полки, К; ТЬ - температура поверхности, К.
Для математического описания второго периода сушки были сделаны следующие допущения:
- содержание влаги в материале в момент начала второго периода сушки равно равновесному влагосодержанию при средней температуре высушенного слоя в конце этапа сублимации (V = ^ед);
- механизмом массопереноса является десорбция влаги.
Баланс влаги в высушенной области мате-
дЖ
риала рассчитывался по формуле -=
дг
= К(1Г - У), где IV- средняя влажность высушенного слоя, кг/кг; К - кинетический коэффициент, с-1; 1Г - равновесное содержание влаги во втором периоде сушки, кг/кг.
Для расчета кинетического коэффициента К
60 Д.,
использовалось уравнение K = ■
'eff
à.
где
Бе^- - эффективный коэффициент диффузии, м2/с; dz - диаметр пор, м.
Для расчета теплового баланса в высушенном слое использовалось уравнение
-Л. ^ + (1 -ф4Л ^ + £ АН ^ ™ = о, * д, у ' дг к М„ дг '
где ке - эффективная теплопроводность высушенного материала, Вт/(м-К); АН - теплота адсорбции, Дж/моль; 5 - общее число компонентов в смеси; с5 - удельная теплоемкость адсорбента, Дж/(кг-К).
Начальные условия для решения уравнений второго периода сушки: Г = ^ед, ( = 0,
Т = Т(х, 0), t = 0.
Граничные условия для решения уравнений второго периода сушки:
Т = Т0, х = 0, t > 0, — = 0, х = Ь, t > 0.
дх
Решение дифференциальных уравнений осуществлялось с помощью метода Мак-Кор-мака.
Алгоритм расчета кинетики вакуумной сублимационной сушки
Для реализации расчета кинетики процесса вакуумной сублимационной сушки описанная
математическая модель была оформлена в программный код в Visual Studia Community на языке программирования C# (собственное ПО). Разработанная программа для расчета кинетики сушки была зарегистрирована как полезная для ЭВМ [10]. На рисунке 2 представлена блок-схема алгоритма расчета кинетики процесса сушки.
С помощью разработанной программы пользователь может рассчитать кинетику вакуумной сублимационной сушки материала - получить значения изменения температуры и влагосодержания материала в процессе сублимационной сушки, предварительно задав физические свойства материала в файле конфигураций. На рисунке 3 представлена последовательность действий, которую должен соблюсти пользователь при работе в программе: 1 - задать физические свойства материала в файле конфигураций c расширением YAML «Имя.уат1.» (данный файл служит для сериа-лизации данных и позволяет хранить сложно-организованные данные в компактном и читаемом виде); 2 - запустить программу для расчета кинетики вакуумной сублимационной сушки и выбрать один из предложенных файлов конфигураций (программа считывает все
Нет
С
Начало
3
Ввод значений в файл конфигураций: s = 0,000001, константы, начальные и граничные условия
Блок рассчитываемых констант
Расчет первого периода сушки Tmi, Wi
Конец
i L.
Экспорт данных в Excel
Да
Нет
Рис. 2. Алгоритм расчета кинетики вакуумной сублимационной сушки Fig. 2. An algorithm for calculating the vacuum freeze-drying kinetics
H HomogeneousDrying yaml ¿3
13 13 20 21 22
23
24
LJU: 0.08 Weg: 0.04 dt: 100 dx: 0.08
Ro sol Ro i ce Çp_sol Cp ice
1040 919. 4 4060 1943 k_sol: 0.546 k_ice: 2.39 rp: 0.00000000421 Rgaz: 8.314 M: 0.018 dz: 0.0098 e_por: 0.667 Eo: 6000 Dso: 0.0Q0000001 Wo: 0.3103
1091.82
0.0017
1111.02
const a: const b: constc: eps: 0.001 P_I: 100 P_II: _00.01868 P_III: 100.06315 P_IV: 100.0412 P_V: 100.02534 T I: 273.15 T_II: 278.15 Till: 283.15 T_IV: 288.15 T_V: 293.15)
C:\Users\4okliova\Desk...
Hello!
Choose corfig-file: 1 - HomogeneousDrying
Выбра ь C:\Users\Mokhova\Desktap\Dr^ing_Mod\Hamageneous_Drying\Mater a Dry n...
Time = 83200; T. Shelf = 293 15 T_m = 293j05658781196836176579514204 Wp_ new = 3,0411956914071042996
Time = 83300; Shelf = 293 15 T_m = 293,057B0696492239786304961921 Wp_ new = 3,9411860781072096238
Time = 83400; T_ Shelf = 293 15 T_m = 293,05901052529121957575418830 Wp_ new = 3,0411765429383991190
Time = 83500; T_ Shelf = 293 15 T_m = 293,06019868736077287398507932 Wp_ new = 3,0411670852498770768
Time = 83600; T. Shelf = 293 15 T_m = 293,06137164307983459414564296 Wp_ new = 3,0411577043965566438
Time = 83700; Shelf = 293 15 T_m = 293,06252958208674376222489312 Wp_ new = 3,0411483997390046436
Time = 83800; T_ Shelf = 293 15 T_m = 293,06367269173585061893513110 Wp_ new = 3,0411391706433871884
Time = 83900; T. Shelf = 293 15 T_m = 293,06480115712368591517933418 Wp_ new = 3,0411300164814158170
Time = 84000; Shelf = 293 15 T_m = 293,06591516111485303071191316 Wp_ new = 3,0411209366302941592
Time = 84100; T_ Shelf = 293 15 T_m = 293,06701488436764545319539763 Wp_ new = 3,0411119304726652864
Time = 84200; T_ Shelf = 293 15 T_m = 293,06810050535939213912568155 Wp_ new = 3,0411029973965596262
Time = 84300; T. Shelf = 293 15 T_m = 293,06917220041153326230462382 Wp_ new = 3,0410941367953434128
Time = 84400; Shelf = 293 15 T_m = 293,07023014371442883968588920 Wp_ new = 3,0410853480676677296
Time = 84500; T_ Shelf = 293 15 T_m = 293,07127450735190270851265853 Wp_ new = 3,0410766306174180782
Time = 84600; T. Shelf = 293 15 T_m = 293,07230546132552431270883731 Wp_ new = 3,0410679838536645678
Time = 84700; Shelf = 293 15 T_m = 293,07332317357863074048313751 Wp_ new = ,041059407190612 5746
Time = 84800; T_ Shelf = 293 15 T_m = 293,07432781002009143906227518 Wp_ new = 3,0410509000475539652
Time = 84900; T. Shelf = 293 15 T_m = 293,07531953454781801638977980 Wp_ new = 3,0410424618488188934
Time = 85000; T. Shelf = 293 15 T_m = 293,07629850907202152351470611 Wp_ new = 3,0410340920237280388
Time = 85100; Shelf = 293 15 T_m = 293,07726489353821959525392424 Wp_ new = 3,0410257900065454376
Time = 85200; T_ Shelf = 293 15 T_m = 293,07821884594999581054658274 Wp_ new = 3,0410175552364317646
Time = 85300; T. Shelf = 293 15 T_m = 293,07916052239151361773363368 Wp_ new = 3,0410093871573981792
Time = 85400; Shelf = 293 15 T_m = 293,08009007704978715379271873 Wp_ new = 3,0410012852182606038
Time = 85500; T_ Shelf = 293 15 T_m = 293,08100766223671127034288348 Wp_ new = 3,0409932488725945662
CalculateDry.
TOTAL CALCULATION TIME: 381ms
EXPORT DATA TO EXCEL: 856
Press any key for exit
Рис. 3. Пример работы программы для расчета кинетических характеристик процесса
вакуумной сублимационной сушки
Fig. 3. An example of the program for calculating the kinetic characteristics of the vacuum freeze-drying process
файлы с расширением УЛМЬ, которые хранятся в папке с проектом); 3 - запустить расчет, при этом значения расчетного времени, температуры полки, а также температуры и влаго-содержания материала выведутся в консоль. Полученные данные по завершении расчета автоматически записываются в Ехсе1-файл, который сохраняется в папке проекта.
Разработанная программа имеет перспективы для применения в области математического моделирования процессов и аппаратов химической технологии, а также в области биотехнологии, так как позволяет рассчитывать и подбирать технологические режимы вакуумной сублимационной сушки.
Построение 3D-геометрии лиофилизатора
Для создания CFD-модели распределения водяных паров по объему рабочей камеры сначала необходимо построить 3D-геометрию сублимационного оборудования. За основу создания 3D-геометрии была принята вакуумная
сублимационная установка ScanVac Coolsafe 100-9 (Дания).
Для упрощения расчетов и уменьшения количества ячеек в расчетной сетке 3D-геометрия вакуумной сублимационной установки была представлена в виде камеры сублиматора с выходным отверстием для водяных паров снизу, которое задавалось как область Pressure-Outlet. Внутри камеры находятся пять нагревательных полок с одной чашкой Петри на каждой полке. Область Pressure-Outlet выступала в качестве упрощенного конденсора с температурой 173 К и давлением 100 Па, так как детальное моделирование змеевика конденсора в данной работе не проводилось. На рисунке 4 представлена 3D-геометрия вакуумной сублимационной камеры ScanVac Coolsafe 100-9, которая создавалась в соответствии с реальными размерами экспериментальной установки ScanVac Coolsafe 100-9 и ее комплектующими. Геометрические размеры рабочей камеры приведены в таблице 1.
Область материала задавалась как mass-flow-inlet, именно к ней прикреплялись UDF-файлы температуры и влагосодержания.
а) б) в)
Рис. 4. 3Б-геометрия сушильной камеры: а) вид спереди, б) изометрическое изображение, в) срез Fig. 4. The drying chamber 3D-geometry: а) front view, б) isometric image, в) slice
CFD-модель распределения водяных паров в сушильной камере
Для решения задачи газодинамики при моделировании потока водяного пара в камере сублиматора необходимо провести процедуру настройки модели. В работе использовалась многофазная модель Эйлера: делается допущение о том, что расчетная область внутри камеры сублиматора состоит из воздуха (фаза 1), представленного неабсолютным вакуумом, и водяного пара (фаза 2), который выделяется из материала по мере процесса. Режим вязкости был выбран как ламинарный. Также для реализации теплопереноса от нагревательных полок к материалу включалось уравнение энергии.
При CFD-моделировании потока водяного пара в вакуумной сублимационной установке и создании расчетной области был принят ряд следующих допущений:
- в начальный момент времени в камере сублиматора присутствует воздух, так как вакуум неабсолютный, а по мере процесса в камеру из материала выделяются водяные пары;
- модель учитывает первый и второй периоды сушки;
- температура стенки сушильной камеры постоянная;
- в камере сублиматора находятся пять полок, которые индексируются следующим образом: верхняя полка - № 5, нижняя полка - № 1;
- рассматривается случай, когда материал с прикрепленными UDF-файлами находится только на пятой и первой полках.
Представим систему балансовых уравнений, основанных на положениях механики сплошных сред. Полученная система дифференциальных уравнений решается численными методами в рамках определенной геометрии и включает уравнения сохранения массы, импульса и энергии соответственно:
d^ + V{pvYl) = V{pw)-Jn
+ V (pw ) = -VP + V (ти ) + pg,
+ E + p)) = V(XVT),
где p - плотность газовой смеси, кг/м3; Y - массовая доля i-го компонента, кг/кг; V - оператор
Таблица 1 Table 1
Геометрические размеры сушильнои камеры и комплектующих
Geometric dimensions of the drying chamber and accessories
Вакуумная камера Нагревательные полки Чашки Петри Материал
h, м 0,43 h, м 0,01 h, м 0,014 h, м 0,01
D^ м 0,277 D, м 0,25 D^ м 0,099 D, м 0,09
DB, м 0,272 H, м 0,045 Du, м 0,093 - -
Набла; v - вектор скорости газа, м/с; J - диффузионный поток i-го компонента; p - давление, Па; ти - тензор вязких напряжений; g -ускорение свободного падения, м/с2; E - полная удельная энергия, Дж/кг; X - коэффициент теплопроводности, Вт/(м К); Т - температура газовой смеси, К.
В программном пакете Ansys Fluent 17.0 дифференциальные уравнения решаются с использованием метода конечных объемов, суть которого заключается в преобразовании дифференциальных уравнений в частных производных к форме алгебраических уравнений.
Результаты математического моделирования
Совмещенное моделирование использовалось для подбора технологического режима сушки. Математическое моделирование кинетики процесса сушки проводилось для 2 масс. % водного альгината натрия (гидрогель). В соответствии с [11] данный гидрогель является основой при получении материалов медицинского назначения. Толщина слоя принималась
равной 0,01 м, начальное содержание влаги в материале - 0,98 кг/кг. В таблице 2 представлены рассчитанные температурные режимы организации вакуумной сублимационной сушки.
Расчет осуществлялся до тех пор, пока общее время сушки (тсушки) не достигало времени режима (трежима).
На рисунке 5 в качестве примера работы разработанного ПО [10] приведены рассчитанные температурные профили материала и изменения влагосодержания для режима 3.
Работа в режиме 1 позволяет интенсифицировать время сушки за счет быстрого поднятия температуры полки до 293,15 К во втором периоде сушки, однако это может негативно сказаться на качестве получаемого продукта. В отличие от режима 1 при работе в режиме 2 поднятие температуры до 293,15 К происходило более плавно. Работа в режиме 3 наиболее продолжительная по времени за счет поэтапного повышения температуры полки во втором периоде сушки. Такие условия обеспечивают постепенный прогрев термолабильного материала и позволяют получить качественный про-
Таблица 2
Температурно-временные параметры
Table 2
Temperature-time parameters
Режим 1
Тполки K 263,15 293,15
трежима, с 0-22 700 22 700-54 600
Результат Остаточное влагосодержание материала 5,3 %; Тматериала = 289,92 К
Режим 2
Тполки K 263,15 273,15 293,15
трежима, с 0-12 300 12 300-26 100 26 100-63 600
Результат Остаточное влагосодержание материала 4,6 %; Тматериала = 292,48 К
Режим 3
Тполки K 273,15 278,15 283,15 293,15
трежима, с 0-28 000 28 000-54 000 54 000-70 000 70 000-80 000
Результат Остаточное влагосодержание материала 4,2 %; Тматериала = 286,6 К
Режим 4
Тполки K 263,15 273,15 283,15 293,15
трежима, с 0-18 000 18 000-36 000 36 000-54 000 54 000-72 000
Результат Остаточное влагосодержание материала 4,3 %; Тматериала = 289,94 К
Режим 5
Тполки K 273,15 288,15
трежима, с 0-18 000 18 000-36 000
Результат Остаточное влагосодержание 10,6 %; Тматериала = 283,21 К
Режим 6
Тполки K 273,15 278,15 293,15
трежима, с 0-18 000 18 000-36 000 36 000-54 000
Результат Остаточное влагосодержание материала 5,4 %; Тматериала = 287,14 К
Рис. 5. Рассчитанные кинетические кривые для режима 3 Fig. 5. The calculated kinetic curves for mode 3
дукт по завершении процесса. Остаточное вла-госодержание при работе в режиме 3 составило 4,2 %, а при работе в режиме 4 - 4,3 %. При работе в режимах 5 и 6 значения остаточного вла-госодержания высокие, поэтому характеристики получаемых продуктов могут не соответствовать требуемому качеству. Анализ результатов серии вычислительных экспериментов (табл. 2) показал, что режим 3 может быть предпочтительнее, так как при переходе к положительным температурам было достигнуто первое критическое влагосодержание (0,15 кг/кг), следовательно, риски «плавления» и растрескивания материала минимальны.
Разработанная CFD-модель распределения водяных паров в объеме вакуумной сублимационной камеры решалась совместно с описанной выше математической моделью кинетики процесса сушки посредством передачи значений массового потока испаряемой влаги и температуры материала через файлы UDF в заданные интервалы времени (через каждые 1 000 с). Данные, записанные в UDF-файлы, использовались в качестве граничных условий для полок 1 и 5.
На рисунке 6 представлен пример организации UDF-файлов для потока испаряемой влаги (UDF-mass_flow_rate) и температуры материала (UDF-temperature_flow).
Представленные на рисунке 6 UDF-файлы были записаны по результатам работы программы для расчета кинетики вакуумной сублимационной сушки для выбранного режима 3.
CFD-модель позволила прогнозировать давление в объеме камеры сублиматора. Полученные значения давления паров над поверхностью материала усреднялись и передавались в кинетическую модель в качестве граничных условий для следующего временного интервала.
Профили распределения давления и скорости потока водяного пара в объеме рабочей камеры в разные моменты времени, рассчитанные для режима 3 (см. http://www.swsys.ru/up-loaded/image/2021 -3/2021 -3 -dop/3jpg), показывают, что максимальное давление (Р = 104,729 Па) на протяжении всего процесса достигается в момент времени 30 000 секунд, что может быть связано с интенсивным выделением водяных паров с поверхности материала в первом периоде сушки. Из полученных результатов также следует, что профиль давления уменьшается сверху вниз по высоте вакуумной камеры. Наибольший градиент давления наблюдается для материала на 5-й полке (АР = 1,119 Па), а наименьшее изменение давления происходит на 1-й полке (АР = 0,342 Па). Чем больше расстояние между конденсором и материалом на полке, тем выше градиент давления. Однако изменения давления незначительны, и для данной геометрии происходит быстрое выравнивание его профиля.
В момент времени 1 000 секунд поток водяных паров незначительный и скорость потока минимальная (у = 0,28249 м/с). Это связано с тем, что данный интервал времени соответ-
UDF-mass_flow_rate
#include "udf.h"
DEFINE_PROFILE(mass_flow_rate,th,i) {
face_t f;
real flow_time = CURRENT_TIME;
begin_f_loop(f,th) {
if(flow_time <= 1000.) F_PROFILE(f,th,i) = 0.0000136; else if(flow_time <=2000) F_PROFILE(f,th,i) = 0.0000273; else if(flow_time <=3000) F_PROFILE(f,th,i) = 0.0000410; else if(flow_time <=4000) F_PROFILE(f,th,i) = 0.0000546; else if(flow_time <=5000) F_PROFILE(f,th,i) = 0.0000683; else if(flow_time <=6000) F_PROFILE(f,th,i) = 0.0000820; else if(flow_time <=7000) F_PROFILE(f,th,i) = 0.0000957; else if(flow_time <=8000) F_PROFILE(f,th,i) = 0.0001093; else if(flow_time <=9000) F_PROFILE(f,th,i) = 0.0001230; else if(flow_time <=10000)
F_PROFILE(f,th,i) = 0.0001367; }
end_f_loop(f,th); }
UDF-temperature_flow
#include "udf.h"
DEFINE_PROFILE(temperature_flow,th,i) {
face_t f;
real flow_time = CURRENT_TIME;
begin_f_loop(f,th) {
if(flow_time <= 1000.) F_PROFILE(f,th,i) = 259.410; else if(flow_time <=2000) F_PROFILE(f,th,i) = 266.167; else if(flow_time <=3000) F_PROFILE(f,th,i) = 270.231; else if(flow_time <=4000) F_PROFILE(f,th,i) = 271.989; else if(flow_time <=5000) F_PROFILE(f,th,i) = 272.728; else if(flow_time <=6000) F_PROFILE(f,th,i) = 273.003; else if(flow_time <=7000) F_PROFILE(f,th,i) = 273.100; else if(flow_time <=8000) F_PROFILE(f,th,i) = 273.133; else if(flow_time <=9000) F_PROFILE(f,th,i) = 273.144; else if(flow_time <=10000)
F_PROFILE(f,th,i) = 273.147; }
end_f_loop(f,th); }
б;
Рис. 6. Пример организации UDF-файлов: а) для потока испаряемой влаги, б) для температуры материала
Fig. 6. An example of the UDF-files organization: a) for the evaporated moisture flow,
б) for the material temperature
ствует периоду прогрева материала. Затем скорость потока водяного пара постепенно растет и достигает максимума (v = 2,802 м/с) в момент времени 30 000 секунд. Скорость водяных паров сильно растет в области Pressure-Outlet с давлением 100 Па.
Значения давлений, полученные в результате CFD-моделирования, фиксировались над 1-й и 5-й полками. Затем эти значения передавались в математическую модель кинетики процесса сушки как граничные условия для уточненного расчета кинетических характеристик отдельно для верхней и нижней полок с целью исследования скорости сушки материалов, находящихся на разном расстоянии от конденсора. Результаты расчета кинетики сушки для нижней и верхней полок для массовых потоков 1410-6, 4110-6, 68 10-6 и 9510-6 кг/с по-
казывают незначительные изменения в давлениях между нижней и верхней полками и ведут к несущественным отличиям в кинетике процесса сушки для данных полок (см. http:// www.swsys.ru/uploaded/image/2021 -3/2021-3-dop/4.jpg). Однако предполагается, что увеличение масштаба лиофилизатора может привести к значительному росту градиента давления и существенно повлиять на однородность и качество получаемого продукта. Данное предположение будет проверено в последующих вычислительных работах, при моделировании крупногабаритного лиофилизатора.
Для значений массового потока 4110-6, 68 10-6 и 95 10-6 кг/с время сушки увеличилось на 850, 1 750 и 2 950 секунд соответственно в сравнении с исходным массовым потоком 14 10-6 кг/с.
Заключение
В статье представлен метод объединения двух математических моделей, используемый для прогнозирования и исследования влияния градиентов давления на скорость процесса сушки материалов, находящихся на разном расстоянии от конденсора. Несмотря на то, что для данной модели вакуумной сублимационной установки (ScanVac Coolsafe 100-9) не отмечено существенных изменений в кинетике сушки материалов на нижней и верхней полках, предполагается, что при масштабировании оборудования градиент давления может резко возрасти. Следовательно, материалы на разных полках будут достигать различной глубины
обезвоживания, что может сказаться на однородности партии.
Математическая модель, разработанная для расчета кинетики вакуумной сублимационной сушки (собственное ПО), позволила рассчитать технологические режимы для 2 масс. % водного альгината натрия. Кроме того, совместное использование данной математической модели и модели, реализованной в Ansys Fluent 17.0, дает возможность прогнозировать распределение водяных паров в рабочей камере, что является важным практическим аспектом, так как сокращаются энергетические затраты и время сушки при реализации технологического процесса.
Литература
1. Wang W., Hu D., Pan Y., Niu L., Chen G. Multiphase transport modeling for freeze-drying of aqueous material frozen with prebuilt porosity. International Journal of Heat and Mass Transfer, 2018, vol. 122, pp. 1353-1365. DOI: 10.1016/j.ijheatmasstransfer.2018.02.054.
2. Диденко А.А., Ершова А.Н., Гордиенко М.Г., Меньшутина Н.В. Моделирование вакуумной сублимационной сушки // Программные продукты и системы. 2011. № 4. С. 170-173.
3. Barresi A.A., Pisano R., Rasetto V., Fissore D., Marchisio D.L. Model-based monitoring and control of industrial freeze-drying processes: effect of batch nonuniformity. Drying Technology, 2010, no. 28, pp. 577-590. DOI: 10.1080/07373931003787934.
4. Barresi A.A., Rasetto V., Marchisio D.L. Use of computational fluid dynamics for improving freeze-dryers design and process understanding. Part 1: Modelling the lyophilisation chamber. European Journal of Pharmaceutics and Biopharmaceutics, 2018, no. 129, pp. 30-44. DOI: 10.1016/j.ejpb.2018.05.008.
5. Aktas M., Sozen A., Amini A., Khanlari A. Experimental analysis and CFd simulation of infrared apricot dryer with heat recovery. Drying Technology, 2017, vol. 35, pp. 766-783. DOI: 10.1080/07373937.2016. 1212871.
6. Ganguly A., Nail S.L., Alexeenko A.A. Rarefied gas dynamics aspects of pharmaceutical freeze-drying. Vacuum, 2012, vol. 86, no. 11, pp. 1739-1747. DOI: 10.1016/j.vacuum.2012.03.025.
7. Zhu T., Moussa E.M., Witting M., Zhou D., Sinha K., Hirth M., Gastens M., Shang S., Nere N., So-mashekar S.G., Alexeenko A., Jameel F. Predictive models of lyophilization process for development, scale-up/tech transfer and manufacturing. European Journal of Pharmaceutics and Biopharmaceutics, 2018, vol. 128, pp. 363-378. DOI: 10.1016/j.ejpb.2018.05.005.
8. Nastaj J.F., Witkiewicz K. Mathematical modeling of the primary and secondary vacuum freeze drying of random solids at microwave heating. International Journal of Heat and Mass Transfer, 2009, vol. 52, pp. 4796-4806. DOI: 10.1016/j.ijheatmasstransfer.2009.06.015.
9. Nastaj J.F., Ambrozek B. Modeling of vacuum desorption in freeze-drying process. Drying Technology, 2005, vol. 23, no. 8, pp. 1693-1709. DOI: 10.1081/DRT-200065117.
10. Мохова Е.К., Пальчикова В.В., Гордиенко М.Г., Меньшутина Н.В. Программа для расчета кинетики вакуумной сублимационной сушки материалов: Свид. о регистр. ПрЭВМ. № 2020665375. Рос. Федерация, 2020.
11. Мохова Е.К., Пальчикова В.В., Соколова Е.А., Гордиенко М.Г. Разработка технологии получения материалов для гемосорбции и культивирования клеток // Успехи в химии и химической технологии. 2020. Т. XXXIV. № 11. С. 7-9.
Software & Systems Received 25.02.21
DOI: 10.15827/0236-235X.135.466-476 2021, vol. 34, no. 3, pp. 466-476
A mathematical model of vacuum freeze-dryer with uneven vapor distribution
over the chamber volume
E.K. Mokhova 1, Graduate Student, lizi5263@yandex.ru
M.G. Gordienko 1, Dr.Sc. (Engineering), Associate Professor, gordienko@muctr.ru 1 D. Mendeleyev University of Chemical Technology of Russia, 125480, Moscow, Russian Federation
Abstract. A vacuum freeze-drying process is used to dehydrate thermolabile substances, and is widely used in biotechnology, pharmaceutical and food industries to obtain dry tissues, drugs and products without losing their structural integrity and biological activity. Nowadays, there is a need to study the influence of pressure gradients in the vacuum chamber on the drying rate of materials located at different distances from the condenser. Uneven removal of solvent vapor from the surface of the material during the drying process can lead to significant differences in drying rates in different parts of the chamber.
Within the framework of this work, the authors jointly used two models: a mathematical model for calculating the kinetics of the vacuum freeze-drying process (native software); Computational Fluid Dynamics (CFD) model for predicting water vapor pressure above a material surface (Ansys Fluent 17.0 software package). The data between the models was transferred by transferring the values of the mass flow of evaporated moisture calculated by the drying kinetics model, through UDF files at specified time intervals. UDF data was used as boundary conditions. The obtained distributions of the water vapor partial pressure were transferred to the kinetic model as boundary conditions for the next cycle of iterations. The joint solution of two mathematical models made it possible to obtain the kinetic characteristics of the drying process for 2 mass. % aqueous sodium alginate taking into account the uneven distribution of water vapor over the vacuum chamber volume. The paper shows that for the ScanVac Coolsafe 100-9 dryer model at the considered mass flows, there is a rapid equalization of the pressure profile along the height of the working chamber, which ensures uniform conditions for material drying. Based on the results of a series of computational experiments, the authors give recommendations on the thermal regime of vacuum freeze-drying of medical materials based on alginate hydrogel.
Keywords: vacuum freeze-drying, heat and mass transfer, drying kinetics, CFD, mathematical modeling.
References
1. Wang W., Hu D., Pan Y., Niu L., Chen G. Multiphase transport modeling for freeze-drying of aqueous material frozen with prebuilt porosity. International Journal of Heat and Mass Transfer, 2018, vol. 122, pp. 1353-1365. DOI: 10.1016/j. ijheatmasstransfer.2018.02.054.
2. Didenko A.A., Yershova A.N., Gordienko M.G., Menshutina N.V. Modeling of vacuum freeze-drying. Software and Systems, 2011, no. 4, pp. 170-173 (in Russ.).
3. Barresi A.A., Pisano R., Rasetto V., Fissore D., Marchisio D.L. Model-based monitoring and control of industrial freeze-drying processes: effect of batch nonuniformity. Drying Technology, 2010, no. 28, pp. 577-590. DOI: 10.1080/ 07373931003787934.
4. Barresi A.A., Rasetto V., Marchisio D.L. Use of computational fluid dynamics for improving freeze-dryers design and process understanding. Part 1: Modelling the lyophilisation chamber. European Journal of Pharmaceutics and Bio-pharmaceutics, 2018, no. 129, pp. 30-44. DOI: 10.1016/j.ejpb.2018.05.008.
5. Aktas M., Sozen A., Amini A., Khanlari A. Experimental analysis and CFD simulation of infrared apricot dryer with heat recovery. Drying Technology, 2017, vol. 35, pp. 766-783. DOI: 10.1080/07373937.2016.1212871.
6. Ganguly A., Nail S.L., Alexeenko A.A. Rarefied gas dynamics aspects of pharmaceutical freeze-drying. Vacuum, 2012, vol. 86, no. 11, pp. 1739-1747. DOI: 10.1016/j.vacuum.2012.03.025.
7. Zhu T., Moussa E.M., Witting M., Zhou D., Sinha K., Hirth M., Gastens M., Shang S., Nere N., Somashekar S.G., Alexeenko A., Jameel F. Predictive models of lyophilization process for development, scale-up/tech transfer and manufacturing. European Journal of Pharmaceutics and Biopharmaceutics, 2018, vol. 128, pp. 363-378. DOI: 10.1016/j.ejpb. 2018.05.005.
8. Nastaj J.F., Witkiewicz K. Mathematical modeling of the primary and secondary vacuum freeze drying of random solids at microwave heating. International Journal of Heat and Mass Transfer, 2009, vol. 52, pp. 4796-4806. DOI: 10.1016/j. ijheatmasstransfer.2009.06.015.
9. Nastaj J.F., Ambrozek B. Modeling of vacuum desorption in freeze-drying process. Drying Technology, 2005, vol. 23, no. 8, pp. 1693-1709. DOI: 10.1081/DRT-200065117.
10. Mokhova E.K., Palchikova V.V., Gordienko M.G., Menshutina N.V. A Program for Calculating the Kinetics of Material Vacuum Freeze Drying. Patent RF, no. 2020665375, 2020.
11. Mokhova E.K., Palchikova V.V., Sokolova E.A., Gordienko M.G. Development of technology for producing materials for hemosorption and cultivation of cells. Advances in Chemistry and Chemical Technology, 2020, vol. XXXIV, no. 11, pp. 7-9 (in Russ.).
Для цитирования
Мохова Е.К., Гордиенко М.Г. Математическая модель вакуумной сублимационной сушки с неравномерным распределением паров по объему камеры // Программные продукты и системы. 2021. Т. 34. № 3. С. 466-476. DOI: 10.15827/0236-235X.135.466-476.
For citation
Mokhova E.K., Gordienko M.G. A mathematical model of vacuum freeze-dryer with uneven vapor distribution over the chamber volume. Software & Systems, 2021, vol. 34, no. 3, pp. 466-476 (in Russ.). DOI: 10.15827/0236-235X. 135.466-476.
