Software & Systems

Вычислительная гидродинамика (computational fluid dynamics, CFD) широко используется для решения различных задач. Среди наиболее известных программных пакетов, работающих с CFD, стоит назвать FLUENT, CFX, FLOWIZARD, PHOENICS, STAR-CD, OPEN FOAM [1, 2].

Вычислительная гидродинамика представляет собой подход, который использует численные методы для решения задач, связанных с описанием движения потоков жидкости. Это осуществляется путем разбиения жидкости на мелкие ячейки для создания расчетной области, называемой сеткой, а затем применения к ним подходящего алгоритма для решения различных уравнений движения [3]. Дифференциальные уравнения Навье–Стокса, содержащие законы сохранения массы, импульса и энергии, являются основой CFD и решаются с использованием таких методов дискретизации, как метод конечных разностей, метод конечных элементов в сочетании с численными методами (например метод Particle-in-Cell (PIC – метод частиц в ячейках) для решения задач гидродинамики с одной, двумя и тремя пространственными переменными), метод Marcer-and-Cell (MAC – метод маркеров и ячеек), произвольные методы Лагранжа–Эйлера и модели k-e-турбулентности и так далее [1, 3]. Решение осуществляется для каждой ячейки, в результате чего формируется поле течения.

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

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

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

Обобщенный подход к моделированию, масштабированию и оптимизации биореакторов

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

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

Сначала анализируются физико-химические параметры культуры клеток/бактерий исследуемого процесса, на основе чего моделируются кинетические зависимости роста клеток, продукта и расхода питательных веществ. Затем с использованием по- лученных данных составляется математическая модель биореактора, учитывающая процессы гидродинамики, тепло- и массопереноса и явления стресса. На следующем этапе рассчитываются параметры работы лабораторного реактора, в том числе определяются скорости движения клеток/бактерий в каждой точке пространства реактора. В качестве подхода для определения адекватности разработанной модели и использующегося для масштабирования реактора был предложен метод PIV (Particle Image Velocimetry), метод цифровой трассерной визуализации, позволяющий измерять векторные поля скорости в выбранном сечении потока газа или жидкости. Одним из достоинств данного метода является возможность измерения в широком динамическом диапазоне скоростей (порядка 500:1), что делает его весьма привлекательным для исследования сложных турбулентных потоков. В зависимости от результата проверки достоверности математического моделирования происходит либо переход к расчету реактора промышленных размеров, либо корректировка уравнений модели и возврат на стадию моделирования.

Частичная реализация данного обобщенного подхода была рассмотрена на примере моделирования, расчета и оптимизации процесса получения молочной кислоты (одной из органических кислот, получаемых биотехнологическим путем) в биореакторе периодического действия при учете воздействия механического стресса на бактерии Lactobacillus casei (далее частицы) с использованием средств вычислительной гидродинамики CFD.

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

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

–      все бактерии (частицы) имеют сферическую форму и один диаметр (принято 1 мкм);

–      вероятность столкновения частиц в объеме аппарата равна 0;

–      все столкновения частиц со стенками аппарата являются абсолютно упругими, то есть потери энергии отсутствуют;

–      в связи с интенсивным перемешиванием достигается равномерное распределение частиц по всему объему биореактора;

–      плотность, вязкость и коэффициент диффузии культуральной среды постоянны и не зависят от времени и координат в биореакторе;

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

Создание геометрии биореактора и расчетной сетки

На первом этапе моделирования была построена геометрическая модель периодического биореактора смешения. Геометрия емкостного биореактора создана в Gambit, программе-генераторе расчетных сеток, включенной в пакет FLUENT 6.1.22. Все размеры соответствуют аналогичным размерам в реальном биореакторе.

Существует несколько разновидностей расчетных сеток в объемах, самые распространенные из них – гексагональная и тетраэдрическая. На данном этапе была выбрана тетраэдрическая разновидность для большей изотропности рабочей сетки (рис. 2).

Основной параметр расчетной сетки при построении – ее плотность. Чем больше ячеек в сетке, тем меньше ошибка в каждой ячейке, а следовательно, и общая ошибка модели. С другой стороны, плотность расчетной сетки влияет на время сходимости решения, и при большом количестве ячеек для сходимости могут потребоваться значительные временные ресурсы, что зачастую неприемлемо при решении сложных задач. Из-за этого приходится прибегать к различным способам уменьшения числа ячеек, основным из которых является снижение плотности сетки в областях, где поведение потоков в достаточной мере прогнозируемо или где можно пренебречь отклонениями.

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

Поскольку биореактор с перемешивающим устройством представляет собой аппарат со сложной структурой потоков, особое внимание при построении расчетной сетки было уделено мешалке: в связи с наличием мелких деталей плотность сетки на ее поверхности была увеличена. В свою очередь, плотность сетки на поверхности аппарата была снижена для обеспечения приемлемого времени сходимости расчета. Это было реализовано путем задания специальной функции по рабочему объему биореактора. Параметры функции – минимальный и максимальный размеры ячеек, скорость увели- чения размеров ячеек. В результате выбрана расчетная сетка, представленная на рисунке 2, где параметры функции размера имеют следующие значения: минимальный размер ячеек – 1 мм, максимальный – 10 мм, скорость роста – 1,1. Число ячеек в объеме биореактора составило около 110 000.

Моделирование биореактора смешения и результаты расчетов

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

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

.                                           (1)

Вращающаяся система отсчета (Rotating Reference Frame) была использована для моделирования вращения мешалки. Мешалка будет иметь нулевую скорость, скорость вращения рабочего объема вокруг мешалки равна постоянной скорости вращения мешалки в реальном аппарате.

Оценка влияния механического стресса в биореакторе проводилась на основе скорости диссипации энергии ε [7]. Выражение для скорости диссипации энергии имеет вид

,    (2)

где ui – реальная скорость;  – средняя скорость; uj – пульсационная скорость.

На основе теории изотропной турбулентности, разработанной А.Н. Колмогоровым, турбулентность объясняется наличием большого количества вихрей. Первичные неустойчивые вихри большого размера, разрушаясь, образуют вихри меньшего размера, при этом имеет место диссипация энергии. Размер вихрей определяет так называемый масштаб турбулентности [7–9]. Для нахождения реакции микроорганизмов на гидродинамический стресс необходимо сравнить их размер с характеристическим размером турбулентных вихрей:

,                                                            (3)

где ν – кинематическая вязкость.

В пакете Fluent для математического описания двухфазных течений чаще всего используют два метода: дискретно-траекторный Эйлер–Лагранжа (Euler–Lagrange) и многожидкостной Эйлер–Эйлера (Euler–Euler). В первом методе рассчитываются траектории и характеристики индивидуальных частиц в определенные промежутки времени, то есть для частиц используется метод Лагранжа, а для жидкой фазы метод Эйлера. Во втором методе дисперсная фаза (частицы) рассматривается как сплошная среда (континуум), а общее течение – как смесь двух или более жидкостей [10].

Для данного процесса объемное содержание дисперсной фазы (бактерий) в относительно малом объеме потока не превышает 10–12 %, внутренний объем биореактора представляет собой сплошную среду. Такая система рассматривается как однофазная, поток является изотропным, поэтому для дисперсной фазы (бактерий) был применен метод Эйлер–Эйлера. Для многофазного потока была применена многофазная модель Эйлера, решающая совместно уравнение сохранения массы:

,                                              (4)

,                            (5)

и уравнения движения:

                      (6)

                             (7)

для каждой из фаз.

Индексы 1 и 2 обозначают сплошную (культуральная среда) и дисперсную (бактерии) фазы соответственно; t – время (с); uq – вектор скоростей q-й фазы (м/с); p – давление в биореакторе (Па); ρq – плотность q-й фазы (кг/м3); g – ускорение свободного падения (м/с2); ps – давление твердой фазы (Па); F – внешние силы (подъемная сила и др.) (Н); Kpq – коэффициент, характеризующий обмен импульсом между фазами p и q, определяется при помощи эмпирических соотношений; αq – объемная доля q-й фазы; τ – тензор напряжений q-й фазы (Па); μq – динамическая вязкость q-й фазы (Па⋅с); λq – объемная вязкость q-й фазы, характеризует деформацию объемного сжатия (Па⋅с).

Для турбулентного режима движения в биореакторе (рассчитанные значения критерия Рейнольдса составили от 3 000 до 15 000 в зависимости от скорости вращения мешалки) во Fluent была выбрана модель k–ε как одна из наиболее распространенных моделей турбулентности, хорошо отображающих потоки, скорость которых на порядок ниже скорости звука. Альтернативная модель k–ω с большей точностью описывает пристеночные эффекты, но сходимость по сравнению с выбранной моделью k–ε хуже. Кроме того, модель k–ε позволяет вычислять скорости диссипации энергии по аппарату, таким образом показывая уровень механического стресса в аппарате.

Выбранная модель k–ε в применении к двухфазной системе может представляться различными уравнениями. При моделировании использовалась дисперсная модель турбулентности, которую целесообразно использовать при небольших объемах дисперсных фаз. Данная модель представляется следующими уравнениями:

 (8)

 (9)

где kq – кинетическая энергия турбулентности q-й фазы, м2/с2; εq – скорость диссипации кинетической энергии q-й фазы, м2/с3; Gkq – производство кинетической энергии q-й фазы, вызванное градиентом скорости;  – вектор средневзвешенных скоростей q-й фазы;  – влияние дисперсной фазы на сплошную фазу q; σk и σε – числа Прандтля соответственно для k и ε; С1ε, С2ε – эмпирические константы, полученные из опытов с воздухом и водой и являющиеся универсальными для большого числа задач гидродинамики и газодинамики. Значения констант: С1ε = 1,44, С2ε = 1,92.

Граничные условия для данных уравнений:

kq(x0; y0; z0) = kq,0  ,                                                 (10)

eq(x0; y0; z0) = eq,0  .                                                 (11)

Результаты и обсуждение

На основании расчетов были получены профили скоростей диссипации энергии по аппарату. Подтверждено, что наибольшие значения скорости диссипации энергии, выделенные красным и желтым цветом, наблюдаются около мешалки (рис. 3). Кроме того, определены профили характеристических размеров турбулентных вихрей по объему биореактора при различных скоростях вращения мешалки (рис. 4). В качестве примера даны графики для скорости перемешивания 1 000 об./мин. Наименьшие значения характеристических размеров вихрей, оказывающие разрушительное действие на клетки, соответствуют области вблизи мешалки.

Было определено, что бактерии подвергаются воздействию механического стресса при размерах турбулентных вихрей, меньших 10 мкм (что соответствует размерам клеток, 0,1–10 мкм). Зона, где наблюдается стрессовое воздействие, называется стрессовой зоной [11]. Увеличение скорости вращения мешалки приводит к возрастанию объема стрессовой зоны в биореакторе (рис. 5). Было введено понятие удельного объема стрессовой зоны, который равен отношению объема стрессовой зоны к рабочему объему биореактора. На основании полученных результатов построена зависимость удельного объема стрессовой зоны биореактора от скорости вращения мешалки (рис. 6). Анализ полученной зависимости показывает, что при достижении примерно 700 об./мин. начинается резкое возрастание удельного объема стрессовой зоны и, соответственно, самого объема этой зоны. Таким образом, будет наблюдаться значительное увеличение воздействия механического стресса на клетки в биореакторе. Превышение данного значения скорости перемешивания приводит к снижению продуктивности по основному продукту – молочной кислоте.

Таким образом, был предложен обобщенный подход к моделированию, масштабированию и оптимизации биореакторов, один из этапов которого заключается в использовании современного метода цифровой трассерной визуализации PIV, позволяющий определять скорости частиц по их изображениям. Частичная реализация данного подхода была рассмотрена на примере моделирования, расчета и оптимизации процесса получения молочной кислоты бактериями Lactobacillus casei в реакторе смешения периодического действия с использованием программного пакета CFD. При расчете учитывалось влияние одного из аспектов механического стресса, перемешивания, на жизнедеятельность и продуктивность бактерий. В качестве характеристики механического стресса выбрана максимальная скорость диссипации энергии. Было введено понятие удельного объема стрессовой зоны, на основании которого определена предельная скорость перемешивания в данном биореакторе смешения, превышение которой ведет к значительному воздействию механического стресса на бактерии. Данный подход может быть использован для перехода к большим реакторам и поиску оптимальных условий работы.

Литература

1.     Sharma C., Malhotra D., Rathore A.S. Review of computational fluid dynamics applications in biotechnology processes. Biotechnol. Prog., 2011. vol. 27, no. 6, pp. 1497–1510.

2.     Hutmacher D.W., Singh H. Computational fluid dynamics for improved bioreactor design and 3D culture. Trends in Biotechn, 2008, vol. 26, no. 4, pp. 166–172.

3.     CFD online. URL: http://www.cfd-online.com (дата обращения: 06.09.2015).

4.     Kaiser S.C., Löffelholz C., Werner S., Eibl D. CFD for characterizing standard and single-use stirred cell culture bioreactors, in: Minin I. (Eds.), Intech, 2011, pp. 97–122.

5.     Johnson С., Natarajan М., Antoniou С. Verification of energy dissipation rate scalability in pilot and production scale bioreactors using computational fluid dynamics. Biotechnol. Progr., 2014, vol. 30, no. 6, pp. 760–764.

6.     Singh H., Hutmacher DW. Bioreactor studies and computational fluid dynamics. Adv Biochem Eng Biotechnol, 2009, vol. 112, pp. 231–249.

7.     Henzler H.-J. Partical stress in bioreactors. Adv. Biochem. Eng. / Biotechnol., 2000, vol. 67, pp. 35–82.

8.     Колмогоров А.Н. Локальная структура турбулентности в несжимаемой жидкости при очень больших числах Рейнольдса // Докл. АН СССР. T. 30. № 4. С. 299–303.

9.     Hu W., Berdugo C., Chalmers J.J. The potential of hydrodynamic damage to animal cells of industrial relevance: current understanding. Cytotechn, 2011, vol. 63, pp. 445–460.

10.  William J. Kelly. Using computational fluid dynamics to characterize and improve bioreactor performance. Biotechnol. Appl. Biochem., 2008, vol. 49, pp. 225–238.

11.  Boudrant J., Menshutina N.V., Skorokhodov A.V., Guse- va E.V., Fick M. Mathematical modelling of cell suspension in high cell density conditions. Application to L-lactic acid fermentation using lactobacillus casei in membrane bioreactor. Process Biochem. 2005, vol. 40, no. 5, pp. 1641–1647.