Теоретические основы химической технологии, 2022, T. 56, № 3, стр. 273-281

Описание работы биохимического реактора диффузионной моделью

Л. Ю. Александрова a*, А. И. Мошинский a

a Санкт-Петербургский химико-фармацевтический университет
Санк-Петербург, Россия

* E-mail: pisce-capricorn@inbox.ru

Поступила в редакцию 04.02.2022
После доработки 08.02.2022
Принята к публикации 15.02.2022

Полный текст (PDF)

Аннотация

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

Ключевые слова: реактор, устойчивость, логистический закон, число Пекле, закон Олли, диффузионная модель, конвективный перенос

ВВЕДЕНИЕ

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

Логистическое уравнение для описания размножения популяций организмов: W(C) = AC(1 − C/K) получило широкое распространение не только в биологических проблемах [1 с. 157, 2 с. 182]. Т.е. оно приобрело, в известном смысле, фундаментальное значение. Биологические популяции ведут себя заметно сложнее типичных физических систем. Тут все реагирует на мелкие возмущения среды обитания и биосистемы постоянно перестраиваются. На данном этапе исследований качественные характеристики подобных систем выдвигаются на передний план. В этой связи предложен ряд коррекций логистической модели, приближающий ее выводы к практическим наблюдениям. Отмеченные уточнения вносят в модель дополнительные параметры, которые хотя и приближают модель к реальности, делают ее анализ более сложным. Здесь мы используем формулу (Олли для популяции с критическим порогом плотности [3]), актуальную, в частности, для экологических проблем

(1)
$W\left( C \right) = AC(1 - {C \mathord{\left/ {\vphantom {C K}} \right. \kern-0em} K})(C - B),~\,\,\,0 < B < K.$

Здесь C — концентрация особей, A, B, K — постоянные, неотрицательные параметры модели. Величина B как раз и определяет порог плотности (концентрацию), ниже которого ростом организмов можно пренебречь.

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

ПРОТОЧНЫЙ БИОХИМИЧЕСКИЙ РЕАКТОР ИДЕАЛЬНОГО ПЕРЕМЕШИВАНИЯ

Уравнение для описания эволюции популяции вытекает из материального баланса и имеет вид

(2)
$V\frac{{dC}}{{d{{\tau }}}} = AC\left( {1 - {C \mathord{\left/ {\vphantom {C K}} \right. \kern-0em} K}} \right)\left( {C - B} \right) - QC,\,\,\,\,B < K,$
где τ – размерное время. Предполагается, что поток жидкости на входе в биохимический реактор не содержит микроорганизмов. Для анализа проблемы целесообразно привести уравнение (2) к безразмерному виду, например, так

(3)
$G = \frac{C}{K},\,\,\,\,t = \frac{{AK{{\tau }}}}{V},\,\,\,\,b = \frac{B}{K},\,\,\,\,q = \frac{Q}{{AK}}.$

Из соотношений (2), (3) следует, что параметр b меняется в пределах: 0 ≤ b < 1. При b = 0 получим закон роста, названный в работе [3] гиперболическим.

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

В переменных (3) уравнение (2) перепишется следующим образом:

(4)
$\frac{{dG}}{{dt}} = G\left( {1 - G} \right)\left( {G - b} \right) - qG.$

Основной интерес для практики представляют стационарные решения уравнения (4), которые будут зависеть от двух параметров b и q. Какое из стационарных решений будет реализовано, зависит от дополнительного (обычно начального) условия для уравнения (4). Для практических целей важно также, чтобы стационарное решение было устойчивым. Решение уравнения с разделяющимися переменными легко находится в аналитическом виде и поэтому здесь подобными вопросами мы заниматься не будем.

Приравняв в (4) dG/dt нулю, находим уравнение для определения стационарных состояний системы:

(5)
$G[(1 - G)(G - b) - q] = 0.$

Очевидно, что (5) имеет тривиальное решение G = 0. Два других решения, имеющих физический смысл, находятся из квадратного уравнения

(6)
${{G}^{2}} - (1{\text{ }} + b)G + q + b = 0,$

при условии положительности детерминанта этого уравнения. Равенство нулю отмеченного детерминанта происходит при:

(7)
$q = {{{{{(1 - b)}}^{2}}} \mathord{\left/ {\vphantom {{{{{(1 - b)}}^{2}}} 4}} \right. \kern-0em} 4},$

(линия в плоскости b , q) и отделяет вариант с тремя вещественными корнями уравнения (5), от случая одного (тривиального) корня. Первый случай реализуется при q < (1 − b)2/4 , второй — при противоположном знаке неравенства. На линии (7) имеет место вырожденный случай кратности корня уравнений (5) и (6).

При положительном детерминанте уравнения (6) (q < (1 − b)2/4) имеем два вещественных положительных корня

(8)
${{G}_{ \pm }} = \frac{{1 + b \pm \sqrt {{{{\left( {1 - b} \right)}}^{2}} - 4q} }}{2}.$

В нашем распоряжении находится параметр Q или q [в безразмерном виде, см. (3)]. Важен вопрос, при каком значении q из реактора выводится максимальное количество особей (массу) в единицу времени. Т.е. когда величина QC (или qG) достигает максимума. В стационарных условиях имеем:

(9)
$P = \max (qG) = \max \left[ {G\left( {1 - G} \right)\left( {G - b} \right)} \right].$

Методами дифференциального исчисления легко устанавливается, что максимум правой части (9) реализуется в физически интересном интервале G ∈ (b, 1), когда

(10)
$G\left( b \right) = \frac{{1 + b + \sqrt {1 - b + {{b}^{2}}} }}{3}.$

Из соотношения (9) вытекает соответствующее значение параметра q:

(11)
$\begin{gathered} q\left( b \right) = \left[ {1 - G\left( b \right)} \right]\left[ {G\left( b \right) - b} \right] = \\ = \frac{1}{9}\left[ {1 - 4b + {{b}^{2}} + \left( {1 + b} \right)\sqrt {1 - b + {{b}^{2}}} } \right]. \\ \end{gathered} $

Максимальное значение производительности реактора равно:

(12)
$\begin{gathered} P\left( b \right) = q\left( b \right)G\left( b \right) = \\ = \frac{1}{9}\left[ {2G\left( b \right)\left( {1 - b + {{b}^{2}}} \right) - b\left( {1 + b} \right)} \right]. \\ \end{gathered} $

Ситуацию иллюстрирует рис. 1, где приведены случаи графического решения уравнения (5). Стационарные решения уравнения (4) соответствуют точкам пересечения линий G(1 − G)(Gb) с прямыми линиями qG.

Рис. 1.

Графики функций: f = qG, 1q = 0.3; 2q = = 0.181; 3q = 0.166; 4q = 0.1, f = G(1 − G)(Gb) — линия 5 при b = 0.15.

Прямая линия 1 пересекается с кривой 5 только при G = 0, т.е. в этом случае имеется единственное тривиальное решение. Заметим, что точка G = 0 принадлежит линии 5, а не имеющая физический смысл часть кривой 5 ниже оси абсцисс на рис. 1 не показана. Линия 2 соответствует вырожденному (структурно неустойчивому) варианту (катастрофе). При этом параметры b и q связаны зависимостью (7). Прямая 3 проходит через максимальное значение функции G(1 − G)(G − ‒ b), т.е. связана с оптимизацией процесса. Значение q определяется формулой (11). Линия 4 − типичная, когда существует три стационарных решения уравнения (5). Она, в отличие от прямой 2, трансверсально пересекает линию G(1 − G)(Gb). Внизу рис. 1 показаны направления движения точки G, когда она отклонилась от стационарного значения. Стрелки (векторное поле) свидетельствуют (выражают), что решение G+ будет устойчивым, тогда как G − неустойчивым. Данный вывод следует из знака производной dG/dt, который определяется дифференциальным уравнением (4). Знак производной будет положительным при G < G < G+ и отрицательным при G < G и G > G+.

Таблица 1.  

Значения функций g0(a) и g1(a) в выбранных точках

a 7 7.5 8 8.5 9 9.5 10 10.5 11
g0(a) 0.467 0.511 0.539 0.562 0.58 0.596 0.61 0.622 0.633
g1(a) 0.756 0.81 0.843 0.856 0.885 0.899 0.911 0.921 0.93

С ростом параметра q (в размерном виде управляющего параметра объемного расхода Q) из рис. 1 видно, что устойчивая стационарная концентрация G+ уменьшается по величине и в определенный момент (линия 2) происходит “катастрофа типа складка” [4]: решение G+ исчезает и в системе остается только одно устойчивое тривиальное решение G = 0. В таком случае говорят, что популяция вымывается из реактора. Отмеченный переход происходит скачком. Величина скачка, как следует из зависимости (8) при q = (1 − b)2/4 равна (1 + b)/2.

При конкретной реализации процесса соответствующего линии 4 на рис. 1, величина G служит пороговым значением. Если начальное условие удовлетворяет неравенству G(0) < G, то интенсивности зарождения популяции микроорганизмов по закону Олли (1) будет недостаточно для ее выживания в реакторе. Со временем придем к стационарной точке G = 0 (популяция вымывается из системы). Если же G(0) > G, то при t → ∞ реализуется стационарное значение концентрации G = G+. Если задать G(0) = G, то ситуация в системе будет неустойчивой и в зависимости от возмущений параметров она придет либо к точке G = 0, либо к точке G = G+.

Отмеченный выше скачок можно описать при помощи уравнения (4) если принять в нем, что параметры b и q (обычно в нашем распоряжении только q) медленно, в определенном масштабе, меняются со временем (квазистационарное приближение). В таком случае система быстро придет к равновесию и далее будет отслеживать какую-нибудь равновесную ветвь. В нашем примере целесообразно использовать ветвь G+(b, q) (8). Уравнение (4), в соответствии со сказанным, преобразуем к виду

(13)
$\begin{gathered} {{\varepsilon }}\frac{{dG}}{{dt}} = G\left( {1 - G} \right)\left( {G - b} \right) - q\left( {{{\varepsilon }}t} \right)G, \\ b = {\text{const}},\,\,\,\,~q(\varepsilon t) = 0.1 + 0.4{\kern 1pt} \text{[}1 - {\text{exp}}( - \varepsilon t)], \\ \end{gathered} $

характерному для описания релаксационных колебаний и пограничных слоев [5, 6], при 0 < ε $ \ll $ 1. Здесь форма функции qt) выбрана в иллюстративных целях и изменен масштаб времени.

Для случая b = 0.15 ситуацию иллюстрирует рис. 2. При отмеченном выборе функции q(t) и ε = 0.01, практически до точки, определяемой уравнением (7) численное решение (кривая 3) практически совпадает с линией G+ (кривая 1). По достижению отмеченной точки численное решение быстро переходит к другому устойчивому состоянию — тривиальному решению G = 0.

Рис. 2.

Графики функций при b = 0.15: 1G+(b, q), 2 – G+(b, q), 3 – решение уравнения (13) для начального условия G(0) = 0.85 .

На рис. 3 представлены характеристики оптимального вывода популяции из биохимического реактора (10)–(12).

Рис. 3.

Функции параметра b, характеризующие максимальную производительность биореактора. 1G (10), 2q (11), 3P (12).

НЕПРОТОЧНЫЙ РЕАКТОР ИДЕАЛЬНОГО ПЕРЕМЕШИВАНИЯ ПРИ УДАЛЕНИИ ПОПУЛЯЦИИ С ПОСТОЯННОЙ ИНТЕНСИВНОСТЬЮ

Рассматриваемый случай аналог задачи размножения, скажем, рыбы в пруду при постоянной квоте отлова. Анализ варианта размножения особей по логистическому закону, проведен, например, в работах [7, 8]. Там отметили опасность оптимизации процесса, которая может привести к вымиранию популяции.

Обозначим интенсивность “вылова” особей из системы (пруда, биохимического реактора и т.п.) R. Поскольку извлекать особи из системы можно лишь тогда, когда они там имеются, следует уточнить формулу для интенсивности так: RH(C). Здесь H(C) = 1 при C > 0 и H(C) = 0 при C ≤ 0 − функция Хевисайда. Впрочем, отрицательное значение R [функцию H(C) исключаем] можно трактовать как добавление в систему особей с постоянной интенсивностью. В случае интенсивности, определяемой зависимостью RH(C) вместо (2) будем рассматривать следующее уравнение:

(14)
$V\frac{{dC}}{{d{{\tau }}}} = AC\left( {1 - {C \mathord{\left/ {\vphantom {C K}} \right. \kern-0em} K}} \right)\left( {C - B} \right) - RH\left( C \right).$

При приведении уравнения (14) к безразмерному виду можно использовать первые три зависимости (3), а вместо параметра q ввести следующий:

(15)
$r = {R \mathord{\left/ {\vphantom {R {(A{{K}^{2}})}}} \right. \kern-0em} {(A{{K}^{2}})}} = {\text{const}}.$

В безразмерных переменных основное уравнение проблемы примет вид

(16)
$\frac{{dG}}{{dt}} = G\left( {1 - G} \right)\left( {G - b} \right) - rH\left( G \right).$

Из уравнения (16) при dG/dt = 0 получаем уравнение для определения стационарных режимов “отлова” популяции

(17)
$\begin{gathered} G(1 - G)(G - b) = r,~\,\, \Rightarrow \,\,~{{G}^{3}} - \left( {1 + b} \right){{G}^{2}} + \\ + \,\,bG + r = 0,~\,\,\,\left( {G > 0} \right). \\ \end{gathered} $

Ясно, что согласно определению функции Хевисайда, стационарным решением уравнения (16) будет также G = 0. Однако правая часть (16) разрывна, поэтому третий (отрицательный) корень (17) будет посторонним. Как известно, кубическое уравнение имеет три корня. При вещественных значениях параметров в (17) реализуется два варианта: 1) три вещественных корня (в частном случае два из них совпадают) и 2) один вещественный и два комплексно сопряженных корня. Один из вещественных корней (17) всегда будет отрицательным. Это вытекает из того, что при G = 0 левая часть (17) меньше правой: (r > 0), а при G → −∞ левая часть неограниченно возрастает, а значит, с определенного значения G, станет больше правой. Непрерывность функций в (17) доказывает сформулированное утверждение.

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

(18)
$3{{G}^{2}} - 2(1{\text{ }} + b)G + b = 0.$

Основной интерес для нас представляет корень квадратного уравнения (18), который определяется формулой (10), поскольку при вычислении максимума выражения (9) использовалась функция, совпадающая с постоянной r в (17). Подставляя выражение (10) в соотношение (17), получим уравнение “критической линии”. Достаточно очевидно, что при этом будем иметь выражение (12). Таким образом две рассмотренные задачи имеют общие характеристики.

График функции r(b) представлен на рис. 3 (линия 3). Два вещественных положительных решения уравнения (17) будут, когда параметры r и b находятся внутри криволинейного клина, образованного отрезками осей координат и линией 3. Выше этой линии положительных корней у уравнения (17) не будет, а на самой линии реализуется вырожденный случай совпадения двух корней.

На рис. 4 представлены некоторые возможные варианты реализации стационарных процессов в биохимическом реакторе. Прямые линии 1,2 на рис. 4 пересекают кривую 5 в двух точках, т.е. определяют два стационарных режима: устойчивый и неустойчивый. Эти режимы отмечены на рис. 4 для случая линии 2. Для значения параметра r, соответствующего линии 4 на этом рисунке положительных стационарных решений нет. В этом случае реализуется только тривиальное решение G = 0. Стрелки внизу рис. 4 имеют тот же смысл, что и для рис. 1.

Рис. 4.

Графики функций: f = r, 1r = 0.03; 2r = = 0.08; 3r = 0.116; 4r = 0.135, f = G(1 − G)(Gb) — линия 5 при b = 0.15.

Оптимальное значение “отлова” в этом случае будет также определяться соотношением (12). При соответствующем выводе в формулу (9) следует вместо qG подставить r. Теперь линии с разным значением параметра r будут параллельными прямыми (см. рис. 4) (в предельном случае это ось G). В предыдущем варианте (см. рис. 1) прямые линии определялись тангенсом угла наклона к оси G, равным q. Максимальное значение “отлова” реализуется для прямой, касающейся вершины горба [максимум функции G(1 − G)(Gb)]. На рис. 4 это линия 3.

Если стратегия оптимального вывода популяции из системы в случае моделирования уравнением (4) отделена определенным интервалом концентрации популяции от “катастрофического” значения, то в варианте уравнения (16) оптимум непосредственно связан с бифуркационным поведением системы. Поэтому в случае уравнения (4) максимальный “отлов” микроорганизмов будет устойчивым в отличие от варианта с уравнением (16), когда небольшое возмущение в системе может привести к полному вымиранию популяции. В работах [7, 8] при логистическом законе размножения популяции, слагаемое qG в уравнении (4) трактуют как обратную связь, позволяющую при сохранении максимальной производительности системы добиться устойчивости выведения из нее, по сравнению с жесткой квотой “отлова” r в уравнении (16). В нашем случае эти выводы сохраняются при трактовке массообменной системе как биохимического реактора.

Отметим также, что небольшое отклонение коэффициента q от оптимального значения (11) приводит не к самоуничтожению системы [как это было при небольшом отклонении (возрастании параметра r) от оптимального плана в модели (16)], а лишь к небольшому уменьшению производительности реактора.

УЧЕТ ПРОДОЛЬНОГО ПЕРЕМЕШИВАНИЯ И КОНВЕКТИВНОГО ПЕРЕНОСА

Будем рассматривать стационарную задачу о диффузионном массопереносе в биохимическом реакторе, записанную в безразмерном виде. Уравнение

(19)
$\frac{1}{{{\text{Pe}}}}\frac{{{{d}^{2}}G}}{{d{{z}^{2}}}} - \frac{{dG}}{{dz}} + W\left( G \right) = 0,$

и граничные условия Данквертса [9]

(20)
${{\left. {{{dG} \mathord{\left/ {\vphantom {{dG} {dz}}} \right. \kern-0em} {dz}}} \right|}_{{z = 0}}} = {\text{Pe}} \times G,~\,\,\,{{\left. {{{dG} \mathord{\left/ {\vphantom {{dG} {dz}}} \right. \kern-0em} {dz}}} \right|}_{{z = 1}}} = 0.$

Такая задача рассматривалась в работе [10] для случая логистического закона роста популяции. Здесь используется закон роста Олли (1), причем переменные и параметры задачи (19),(20) приведены к безразмерному виду следующим образом:

(21)
$\begin{gathered} z = {x \mathord{\left/ {\vphantom {x L}} \right. \kern-0em} L},~\,\,G = {C \mathord{\left/ {\vphantom {C K}} \right. \kern-0em} K},~ \\ {\text{Pe}} = {{UL} \mathord{\left/ {\vphantom {{UL} D}} \right. \kern-0em} D},~\,\,a = {{AKL} \mathord{\left/ {\vphantom {{AKL} {\left( {UV} \right)}}} \right. \kern-0em} {\left( {UV} \right)}} = {1 \mathord{\left/ {\vphantom {1 q}} \right. \kern-0em} q}, \\ \end{gathered} $
где D − коэффициент дисперсии, G − безразмерная концентрация микроорганизмов в реакторе; L − длина реактора, Pe − число Пекле; U − скорость переноса биомассы, x и z − размерная и безразмерная координаты вдоль оси реактора соответственно.

Условия (20) относится к частному случаю распространенных граничных условий Данквертса, применяемых, как правило, для описания диффузионного продольного перемешивания веществ в реакторе. Первое выражение (20) определяет условие, при котором отсутствует приток массы (в рассматриваемом случае микроорганизмов) в реактор, а второе указывает на равенство потоков на выходе из реактора, когда за его пределами коэффициент дисперсии равен нулю (D = 0) [9].

БОЛЬШИЕ ЗНАЧЕНИЯ ЧИСЛА ПЕКЛЕ

Рассмотренный выше режим идеального перемешивания является предельным вариантом задачи (19), (20) (при учете нестационарного слагаемого), когда Pe → 0 [9, 11]. Противоположный предельный вариант Pe → ∞ называют режимом идеального вытеснения. Выполнив этот предельный переход в соотношениях (19) и (20), получаем, приняв для функции W закон Олли (1):

(22)
$\frac{{dG}}{{dz}} = aG\left( {1 - G} \right)\left( {G - b} \right),\,\,\,\,G\left( 0 \right) = 0.$

Отметим, что с точностью до обозначения независимой переменной, уравнение (22) совпадает с (4) при q = 0 (taz). Задача (22) имеет единственное тривиальное решение G(z) = 0.

Следует заметить, что предельный переход Pe → ∞ понижает порядок уравнения (19), что является признаком сингулярности соответствующего разложения [5, 6]. Поэтому второе граничное условие (20) опущено. Однако в окрестности точки z = 1 при Pe → ∞ появляется пограничный слой, согласующий решения “внешних” задач типа (22) с нетривиальным решением приближения пограничного слоя при граничном условии dG/dzz=1 = 0 [12]. При тривиальном решении “внешней” задачи, решение “внутренней” задачи (приближение пограничного слоя) также будет тривиальным.

ПРОМЕЖУТОЧНЫЕ ЗНАЧЕНИЯ ЧИСЛА ПЕКЛЕ

Устойчивость работы проточных химических реакторов на основе диффузионной модели рассматривалась в ряде работ (напр. [1114]). Этому фактору уделяется серьезное внимание в литературе. Из предыдущего анализа следует, что при фиксированных значениях параметров a и b, с ростом числа Пекле нетривиальное устойчивое решение (которое только и имеет практическое значение) может исчезнуть [в тех случаях, когда при Pe → 0 оно имелось (см. рис. 1)], т.е. естественно предположить, существование граничного значения ${\text{P}}{{{\text{e}}}_{*}}$(a, b) такого, что при Pe > ${\text{P}}{{{\text{e}}}_{*}}$(a, b) задача имеет только тривиальное решение.

Для проверки вышесказанных положений и определения параметров, характеризующих процесс, были проведены расчеты нелинейной задачи (19), (20) для функции W(G) = aG(1 − G)(Gb). Использовался метод пристрелки [15]. Для рассматриваемой задачи “стрельбу” целесообразно вести от точки z = 1 до точки z = 0, поскольку при “прямом” направлении интегрирования (от точки z = 0 до точки z = 1) ошибки округления нарастают быстрее, что отмечалось, в частности в работе [16 с. 295]. Для удобства введем новую независимую переменную y = 1 − z и сведем уравнение второго порядка (19) к системе обыкновенных дифференциальных уравнений первого порядка. Пусть T = dG/dy, тогда задача (19), (20) принимает вид

(23)
$\begin{gathered} \left\{ \begin{gathered} {{dT} \mathord{\left/ {\vphantom {{dT} {dy}}} \right. \kern-0em} {dy}} = - {\text{Pe}}\left[ {aG\left( {1 - G} \right)\left( {G - b} \right) + T} \right] \hfill \\ {{dG} \mathord{\left/ {\vphantom {{dG} {dy}}} \right. \kern-0em} {dy}} = T \hfill \\ \end{gathered} \right., \\ \left\{ \begin{gathered} {{\left. T \right|}_{{y = 0}}} = 0 \hfill \\ T + {{\left. {{\text{Pe}}G} \right|}_{{y = 1}}} = {\text{0}} \hfill \\ \end{gathered} \right.. \\ \end{gathered} $

При реализации метода пристрелки естественно и целесообразно в качестве параметра стрельбы [недостающего (пробного) граничного условия при z = 1] использовать концентрацию микроорганизмов в выходном сечении биохимического реактора (z = 1): G(1) = Ck. Часто именно эта величина представляет основной интерес для практики.

На рис. 5 представлены типичные кривые, демонстрирующие зависимость концентрации на выходе из реактора Ck , как функцию числа Пекле. В случае q > (1 − b)2/4 при любых значениях Pe реализуется только тривиальное решение. Заметим также, что потеря нетривиальных решений происходит скачком.

Рис. 5.

Концентрация на выходе из реактора как функция числа Pe. 1a = 10, b = 0.25; 2a = 8, b = = 0.15; 3a = 10, b = 0.15.

Характер изменения концентрации микроорганизмов вдоль реактора, описываемый нетривиальным решением задачи (19), (20) при W(G) = = aG(1 − G)(Gb) или задачи (23) показан для различных параметров на рис. 6.

Рис. 6.

Распределение концентрации популяции вдоль реактора. 1a = 10, b = 0.15, Pe = 0.5; 2a = 10, b = 0.15, Pe = 1; 3a = 10, b = 0.25, Pe = 1; 4a = 7, b = 0.15, Pe = 1; 5a = 8, b = 0.15, Pe = 1.5.

АППРОКСИМАЦИЯ (ОПИСАНИЕ) НЕКОТОРЫХ ЛИНИЙ А ПРОСТОЙ АНАЛИТИЧЕСКОЙ ЗАВИСИМОСТЬЮ

Одной из проблем (недостатком) решений задач, полученных сложными, в частности численными алгоритмами, для использования в приложениях является нахождение решения при не приведенных в соответствующих источниках значений параметров. Для инженерной практики обычно привлекательны простые расчетные формулы. Особенно это полезно, когда необходимо использовать функции двух и более переменных. В целях удобной аппроксимации экспериментальных данных в работах [17, 18] предложен так называемый метод асимптотических координат. Суть метода [17, 18] сводится к построению ряда функций одной переменной для аппроксимации двумерной поверхности путем выбора рациональных для этой цели координат и преобразований. Как правило, метод не является точным, а основан на схожем поведении серий кривых типа приведенных у нас на рис. 6. Ясно, что идеи работ [17, 18] могут быть перенесены на аппроксимацию двумерных поверхностей, полученных расчетным (теоретическим) путем. Это полезно, когда расчетные формулы труднообозримы или расчеты достаточно трудоемки.

В иллюстративных целях покажем, как можно получить простую аналитическую формулу для аппроксимации функции G(z, a, b, Pe) в некотором диапазоне независимых параметров. Для конкретности двумерную поверхность будем строить при b = 0.15, Pe = 1. Несколько модифицируя метод [17, 18], построим две функции G(0, a, 0.15, 1) = = g0(a) и G(1, a, 0.15, 1) = g1(a), для чего вычислим по алгоритму (23) значения приведенных функций в ряде точек. Результаты представлены в табл. 1.

Аппроксимируем в рассматриваемом диапазоне переменной a ∈ (7, 11) функции g0 и g1 выражениями

(24)
$\begin{gathered} g0(a) = 0.104{{(a - 6.8)}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}} + 0.424;~ \hfill \\ g1(a) = 0.108{\text{ }}{{(a - 6.8)}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}} + 0.719. \hfill \\ \end{gathered} $

Визуально о точности аппроксимации табличных значений функций g0 и g1 можно судить по данным рис. 7.

Рис. 7.

Функции: 1g0, 2g1, × − табличные данные.

Теперь, следуя [17, 18] построим вспомогательную функцию

(25)
$F\left( {z,a} \right) = \frac{{G(z,a,0.15,1) - g0(a)}}{{g1(a) - g0(a)}}.$

Сущность методики [17, 18] заключается в том, что функции типа F(z, a) во многих случаях с практически приемлемой точностью не зависят от переменной a и могут быть легко аппроксимированы функциями единственной переменной z.

На рис. 8 представлены некоторые кривые, выражающие функцию F(z, a), а также функция одной переменной g(z), близкая к представленным линиям. Функцию g(z) представляет выражение

(26)
$F(z,a) = g(z) = 1.98\ln (1 + z) - 0.375{{z}^{2}}.$
Рис. 8.

Отдельные кривые функции F(z, a): 1a = 7, 2 − a = 9, 3a = 11. Пунктирная линия − функция g(z).

Не достаточно точная аппроксимация, представленная второй функцией g1 (24), приводит к тому, что не выполняется условие F(1, a) = 1 при любых значениях a, вытекающая из определения этой функции (25). Несколько лучше обстоит дело с аппроксимацией, выраженной функцией g0, точнее условие F(0, a) = 0 можно считать выполненным. Ясно, что, улучшив аппроксимацию (например, путем включения в зависимости (24) дополнительных параметров) можно добиться более близкого друг к другу расположения кривых на рис. 8.

Фактически мы заменили серию близких кривых на рис. 8 одной кривой в соответствующих координатах. Она определяется формулой (26), представляющей собой функцию одной переменной. Мы по существу свели описание двумерной поверхности решения задачи в области z ∈ (0, 1), a ∈ (7, 11) тремя функциями одной (24), (26). В результате из зависимости (25) находим

$G\left( {z,a,{\text{ }}0.15,{\text{ }}1} \right) \approx g\left( z \right)[g1\left( a \right) - g0\left( a \right)] + g0\left( a \right).$

Успех методики, предлагаемой здесь и в [17, 18] зависит от такого выбора вспомогательных функций и переменных, что функция типа F(z, t) (25) с приемлемой точностью зависит только от одной переменной, в идеале, эта зависимость точная (в работах [17, 18] приведен пример такой реализации метода).

ЗАКЛЮЧЕНИЕ

Рассмотренные варианты моделирования биохимических реакторов на основе диффузионной модели и ее предельных случаев показали важность анализа качественной перестройки решений соответствующих уравнений, которые встречаются и в экологических проблемах. Исчерпывающему исследованию рассмотренных вопросов, имеющих также значение для экологии, мешает большое количество параметров, влияющих на качественное поведение системы. Например задача (23) содержит три безразмерных параметра a, b, и Pe.

ОБОЗНАЧЕНИЯ

A параметр в законе Олли, м6/(кг c)
a безразмерный параметр
B параметр в законе Олли, кг/м3
b безразмерный параметр, определен в (3)
C концентрация особей, кг/м3
D коэффициент дисперсии, м2/c
F функция, определена в (25)
G безразмерная концентрация особей
H(C) функция Хевисайда
K параметр в законе Олли, кг/м3
L длина реактора, м
P безразмерная максимальная производительность реактора
Q объемный расход через реактор, м3/c
q безразмерный объемный расход, определен в (3)
R массовый расход особей, кг/c
r безразмерный массовый расход особей T = dG/dy
t безразмерное время, определено в (3)
V объем реактора, м3
W интенсивность производства особей, кг/c
x координата в диффузионной модели, м
y = 1 – z
z безразмерная координата в диффузионной модели
ε безразмерный малый параметр
τ время, c
Pe число Пекле
± два значения величины G определено в (8)
* граничное значение числа Пекле

Список литературы

  1. Гордин В.А. Дифференциальные и разностные уравнения. Какие явления они описывают и как их решать. Учебное пособие. М.: Изд. Дом Высшей школы экономики, 2016.

  2. Трубецков Д.И. Введение в синергетику. Хаос и структуры. Изд. 2-е испр. и доп. М.: Едиториал УРСС, 2004.

  3. Свирежев Ю.М. Нелинейные волны, диссипативные структуры и катастрофы в экологии. М.: Наука, 1987.

  4. Алексеев Ю.К., Сухоруков А.П. Введение в теорию катастроф: учебное пособие. Изд. 2‑е, доп. М.: Книжный дом “Либроком”, 2009.

  5. Коул Дж. Методы возмущений в прикладной математике. М.: Мир, 1972.

  6. Найфэ А. Методы возмущений. М.: Мир, 1976.

  7. Арнольд В.И. Обыкновенные дифференциальные уравнения. Издание 3-е, перераб. и доп. Москва: Наука, 1984.

  8. Арнольд В.И. “Жесткие” и “мягкие” математические модели. Изд. 4-ое, стереотипное. М.: Изд-во МЦНМО, 2013.

  9. Мошинский А.И. Моделирование тепломассообменных процессов на основе обобщенных диффузионных уравнений. М.: Изд-во КНОРУС, 2019.

  10. Кафаров В.В., Винаров А.Ю., Гордеев Л.С. Моделирование биохимических реакторов. М.: Лесная промышленность, 1979.

  11. Мошинский А.И. Описание стационарной работы биохимического реактора диффузионной моделью // Инж.-физ. журн. 2017. Т. 90. № 4. С. 814.

  12. Мошинский А.И. О влиянии релаксационных процессов на устойчивость стационарных режимов работы изотермического проточного химического реактора // Инж.-физ. журн. 1991. Т. 61. № 1. С. 129.

  13. Берман В.С., Курдюмов В.Н., Рязанцев Ю.С. Об устойчивости стационарных режимов изотермического проточного химического реактора // Известия АН СССР, серия Механика жидкости и газа. 1985. № 2. С. 179.

  14. Бучин В.А., Ларин О.Б. Стабилизация неустойчивых режимов работы химического реактора вытеснения // Докл. АН СССР. 1983. Т. 271. № 6. С. 1440.

  15. Современные численные методы решения обыкновенных дифференциальных уравнений / Под ред. Холл Дж. и Уатт. Дж. М.: Мир. 1979.

  16. Арис Р. Анализ процессов в химических реакторах. Л.: Химия. 1967.

  17. Полянин А.Д., Альварес-Суарес В.А., Дильман В.В., Рязанцев Ю.С. Метод асимптотической интерполяции экспериментальных данных // Теорет. основы хим. технологии. 1986. Т. 20. № 5. С. 584.

  18. Дильман В.В., Полянин А.Д. Методы модельных уравнений и аналогий в химической технологии. М.: Химия, 1988.

Дополнительные материалы отсутствуют.