Software & Systems

Разработка математических моделей является одной из актуальных задач в области исследования процесса вакуумной сублимационной сушки. Модели строятся на описании тепло- и массообмена процесса сушки, выраженного с помощью дифференциальных уравнений [1, 2]. В результате решения системы дифференциальных уравнений могут быть получены кинетические характеристики процесса сушки – изменение температуры и влагосодержания материала. Такие математические модели позволяют, варьируя рабочие условия процесса (температуру нагревания полок и глубину вакуума в камере сублиматора), подобрать режимы ведения процесса, обеспечивающие снижение энергетических затрат и времени сушки [3].

За последние несколько лет появился большой интерес к моделированию промышленного и лабораторного оборудования с использованием вычислительной гидродинамики (Computational Fluid Dynamics, CFD). Это связано с тем, что CFD-моделирование можно использовать для улучшения конструкций, оптимизации геометрии и оценки влияния различных вариантов устройства оборудования на ход процесса [4–6].

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

Цель данного исследования – разработка и программная реализация математической модели процесса вакуумной сублимационной сушки с учетом неравномерного распределения водяных паров по объему камеры и их влияния на кинетику процесса. Разработанная математическая модель позволяет подбирать технологические режимы ведения процесса для различных материалов, что существенно сокращает время сушки и энергетические затраты.

Объединение математических подходов

Для реализации цели разработанная математическая модель кинетики процесса сублимационной сушки решалась совместно с CFD-моделью распределения водяных паров по объему вакуумной камеры. Данный подход позволил учесть влияние распределения паров растворителя над поверхностью материала на разном удалении от входа в десублиматор на скорость процесса.

Для этого значения изменения температуры и влагосодержания материала, полученные по кинетической модели, записывались в UDF-файлы temperature_flow и mass_flow_rate, которые прикреплялись к соответствующим областям созданной 3D-геометрии сублимационной камеры в качестве граничных условий в среде моделирования Ansys Fluent 17.0. На рисунке 1 представлен алгоритм связи двух моделей.

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

Математическое моделирование кинетики вакуумной сублимационной сушки

При разработке данной математической модели за основу были приняты статьи [8, 9]. Математическая модель является одномерной, то есть перенос тепла и влаги происходит только вдоль оси Х и строится на уравнениях тепло- и массопереноса. В структуру настоящей математической модели входят описания первого и второго периодов сушки.

Для математического описания первого периода сушки были приняты следующие допущения:

-     существует фронт сублимации Xp, который делит образец на замороженную (I) и высушенную (II) области;

-     граница раздела (фронт сублимации) в процессе сушки равномерно двигается вниз, пока вся замороженная свободная влага не удалится из материала толщиной L;

-     на границе раздела фаз выполняется ус­ловие: массовый поток через границу в направлении снизу вверх равен массовому потоку через границу в направлении сверху вниз;

-     высушенная и замороженная области рассматриваются как гомогенные с точки зрения постоянства теплопроводности, плотности и теплоемкости;

-     теплопередача к материалу осуществляется кондуктивно от полки снизу.

Далее представлено математическое описание первого периода сушки.

Для замороженной области (I) уравнение теплопроводности выглядит следующим образом:  где T – температура, К; t – время, с; aeI – температуропроводность, м2/с; X – декартова координата, м; QvI – объемная мощность источников теплоты, Вт/м3; ρeI – плотность, кг/м3; CpeI – удельная теплоемкость, Дж/(кг‧К).

Уравнение теплопроводности для высушенной области (II) следующее:

где ρbuII – насыпная плотность сухого материала, кг/м3; Nw – поток водяного пара в области II, кг/(м2‧с).

Приращение координаты фронта сублима-ции ∆X*n+1 и влагосодержания ∆W*n+1 в безраз-мерном виде рассчитывалось для областей (I) и (II) соответственно:

где q* – безразмерный тепловой поток, Вт/м2; Stev – модифицированное число Стефана; ∆t* – безразмерный шаг по времени; n – позиция во времени.

Расчет смещения фронта сублимации X*n+1 и изменения влагосодержания W*n+1 проводился также в безразмерном виде соответственно:

Перевод влагосодержания в размерный вид осуществлялся по формуле

где W* – безразмерное влагосодержание; Wp – начальное влагосодержание материала, кг/кг; Weq – равновесная влажность в первый период сушки, кг/кг.

Начальное условие для замороженной и высушенной областей материала: TI = TII, 0 ≤ x ≤ L, t = 0.

Представим граничные условия для нижней (I) и верхней (II) частей материала:

T0 = TI, x = 0, t ˃ 0,

TII = TL, x = L, t ˃ 0,

x = X(t), t ˃ 0,

где T0 – температура полки, К; TL – температура поверхности, К.

Для математического описания второго периода сушки были сделаны следующие допущения:

-     содержание влаги в материале в момент начала второго периода сушки равно равновесному влагосодержанию при средней температуре высушенного слоя в конце этапа сублимации (W = Weq);

-     механизмом массопереноса является десорбция влаги.

Баланс влаги в высушенной области материала рассчитывался по формуле  = K(Wr – W), где W – средняя влажность высушенного слоя, кг/кг; K – кинетический коэффициент, с-1; Wr – равновесное содержание влаги во втором периоде сушки, кг/кг.

Для расчета кинетического коэффициента K использовалось уравнение  где Deff – эффективный коэффициент диффузии, м2/с; dz – диаметр пор, м.

Для расчета теплового баланса в высушенном слое использовалось уравнение

где ke – эффективная теплопроводность высушенного материала, Вт/(м‧К); ∆H – теплота адсорбции, Дж/моль; s – общее число компонентов в смеси; cs – удельная теплоемкость адсорбента, Дж/(кг‧К).

Начальные условия для решения уравнений второго периода сушки:

W = Weq, t = 0,

T = T(x, 0), t = 0.

Граничные условия для решения уравнений второго периода сушки:

T = T0, x = 0, t ˃ 0,  = 0, x = L, t ˃ 0.

Решение дифференциальных уравнений осуществлялось с помощью метода Мак-Кормака.

Алгоритм расчета кинетики вакуумной сублимационной сушки

Для реализации расчета кинетики процесса вакуумной сублимационной сушки описанная математическая модель была оформлена в программный код в Visual Studia Community на языке программирования C# (собственное ПО). Разработанная программа для расчета кинетики сушки была зарегистрирована как полезная для ЭВМ [10]. На рисунке 2 представлена блок-схема алгоритма расчета кинетики процесса сушки.

С помощью разработанной программы пользователь может рассчитать кинетику вакуумной сублимационной сушки материала – получить значения изменения температуры и влагосодержания материала в процессе сублимационной сушки, предварительно задав физические свойства материала в файле конфигураций. На рисунке 3 представлена последовательность действий, которую должен соблюсти пользователь при работе в программе: 1 – задать физические свойства материала в файле конфигураций c расширением YAML «Имя.yaml.» (данный файл служит для сериализации данных и позволяет хранить сложноорганизованные данные в компактном и читаемом виде); 2 – запустить программу для расчета кинетики вакуумной сублимационной сушки и выбрать один из предложенных файлов конфигураций (программа считывает все файлы с расширением YAML, которые хранятся в папке с проектом); 3 – запустить расчет, при этом значения расчетного времени, тем-пературы полки, а также температуры и влагосодержания материала выведутся в консоль. Полученные данные по завершении расчета автоматически записываются в Excel-файл, который сохраняется в папке проекта.

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

Построение 3D-геометрии лиофилизатора Для создания CFD-модели распределения водяных паров по объему рабочей камеры сначала необходимо построить 3D-геометрию сублимационного оборудования. За основу создания 3D-геометрии была принята вакуумная сублимационная установка ScanVac Coolsafe 100-9 (Дания).

Для упрощения расчетов и уменьшения количества ячеек в расчетной сетке 3D-геометрия вакуумной сублимационной установки была представлена в виде камеры сублиматора с выходным отверстием для водяных паров снизу, которое задавалось как область Pressure-Outlet. Внутри камеры находятся пять нагревательных полок с одной чашкой Петри на каждой полке. Область Pressure-Outlet выступала в качестве упрощенного конденсора с температурой 173 К и давлением 100 Па, так как детальное моделирование змеевика конденсора в данной работе не проводилось. На рисунке 4 представлена 3D-геометрия вакуумной сублимационной камеры ScanVac Coolsafe 100-9, которая создавалась в соответствии с реальными размерами экспериментальной установки ScanVac Cool­safe 100-9 и ее комплектующими. Геометрические размеры рабочей камеры приведены в таблице 1.

Область материала задавалась как mass-flow-inlet, именно к ней прикреплялись UDF-файлы температуры и влагосодержания.

CFD-модель распределения водяных паров в сушильной камере

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

При CFD-моделировании потока водяного пара в вакуумной сублимационной установке и создании расчетной области был принят ряд следующих допущений:

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

-     модель учитывает первый и второй периоды сушки;

-     температура стенки сушильной камеры постоянная;

-     в камере сублиматора находятся пять полок, которые индексируются следующим образом: верхняя полка – № 5, нижняя полка – № 1;

-     рассматривается случай, когда материал с прикрепленными UDF-файлами находится только на пятой и первой полках.

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

где ρ – плотность газовой смеси, кг/м3; Y – мас-совая доля i-го компонента, кг/кг; ∇ – оператор Набла; – вектор скорости газа, м/с; J – диффузионный поток i-го компонента; p – давление, Па; τkl – тензор вязких напряжений;  – ускорение свободного падения, м/с2; E – полная удельная энергия, Дж/кг; λ – коэффициент теплопроводности, Вт/(м·К); Т – температура газовой смеси, К.

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

Результаты математического моделирования

Совмещенное моделирование использовалось для подбора технологического режима сушки. Математическое моделирование кинетики процесса сушки проводилось для 2 масс. % водного альгината натрия (гидрогель). В соответствии с [11] данный гидрогель является основой при получении материалов медицин-ского назначения. Толщина слоя принималась равной 0,01 м, начальное содержание влаги в материале – 0,98 кг/кг. В таблице 2 представлены рассчитанные температурные режимы организации вакуумной сублимационной сушки.

Расчет осуществлялся до тех пор, пока общее время сушки (τсушки) не достигало времени режима (τрежима).

На рисунке 5 в качестве примера работы разработанного ПО [10] приведены рассчитанные температурные профили материала и изменения влагосодержания для режима 3.

Работа в режиме 1 позволяет интенсифицировать время сушки за счет быстрого поднятия температуры полки до 293,15 K во втором периоде сушки, однако это может негативно сказаться на качестве получаемого продукта. В отличие от режима 1 при работе в режиме 2 поднятие температуры до 293,15 K происходило более плавно. Работа в режиме 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/uploaded/image/2021-3/2021-3-dop/3.jpg), показывают, что максимальное давление (P = 104,729 Па) на протяжении всего процесса достигается в момент времени 30 000 секунд, что может быть связано с интенсивным выделением водяных паров с поверхности материала в первом периоде сушки. Из полученных результатов также следует, что профиль давления уменьшается сверху вниз по высоте вакуумной камеры. Наибольший градиент давления наблюдается для материала на 5-й полке (∆P = 1,119 Па), а наименьшее изменение давления происходит на 1-й полке (∆P = 0,342 Па). Чем больше расстояние между конденсором и материалом на полке, тем выше градиент давления. Однако изменения давления незначительны, и для данной геометрии происходит быстрое выравнивание его профиля.

В момент времени 1 000 секунд поток водяных паров незначительный и скорость потока минимальная (v = 0,28249 м/с). Это связано с тем, что данный интервал времени соответ-ствует периоду прогрева материала. Затем скорость потока водяного пара постепенно растет и достигает максимума (v = 2,802 м/с) в момент времени 30 000 секунд. Скорость водяных паров сильно растет в области Pressure-Outlet с давлением 100 Па.

Значения давлений, полученные в результате CFD-моделирования, фиксировались над 1-й и 5-й полками. Затем эти значения передавались в математическую модель кинетики процесса сушки как граничные условия для уточненного расчета кинетических характеристик отдельно для верхней и нижней полок с целью исследования скорости сушки материа-лов, находящихся на разном расстоянии от конденсора. Результаты расчета кинетики сушки для нижней и верхней полок для массовых потоков 14·10-6, 41·10-6, 68·10-6 и 95·10-6 кг/с по-казывают незначительные изменения в давлениях между нижней и верхней полками и ведут к несущественным отличиям в кинетике процесса сушки для данных полок (см. http://www.swsys.ru/uploaded/image/2021-3/2021-3-dop/4.jpg). Однако предполагается, что увеличение масштаба лиофилизатора может привести к значительному росту градиента давления и существенно повлиять на однородность и качество получаемого продукта. Данное предположение будет проверено в последующих вычислительных работах, при моделировании крупногабаритного лиофилизатора.

Для значений массового потока 41·10-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., 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. Мохова Е.К., Пальчикова В.В., Гордиенко М.Г., Меньшутина Н.В. Программа для расчета кинетики вакуумной сублимационной сушки материалов: Свид. о регистр. ПрЭВМ. № 2020665375. Рос. Федерация, 2020.

11. Мохова Е.К., Пальчикова В.В., Соколова Е.А., Гордиенко М.Г. Разработка технологии получения материалов для гемосорбции и культивирования клеток // Успехи в химии и химической технологии. 2020. Т. XXXIV. № 11. С. 7–9.

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 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., 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.).