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

ВВЕДЕНИЕ

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

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

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

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

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

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

где τ – размерное время. Предполагается, что поток жидкости на входе в биохимический реактор не содержит микроорганизмов. Для анализа проблемы целесообразно привести уравнение (2) к безразмерному виду, например, так

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Прямая линия 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)(G − b). Внизу рис. 1 показаны направления движения точки G, когда она отклонилась от стационарного значения. Стрелки (векторное поле) свидетельствуют (выражают), что решение G₊ будет устойчивым, тогда как G₋ − неустойчивым. Данный вывод следует из знака производной dG/dt, который определяется дифференциальным уравнением (4). Знак производной будет положительным при G₋ < G < G₊ и отрицательным при G < G₋ и G > G₊.

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

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

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

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

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

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

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

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

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

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

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

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

Основной интерес для нас представляет корень квадратного уравнения (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.

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

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

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

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

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

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

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

где D − коэффициент дисперсии, G − безразмерная концентрация микроорганизмов в реакторе; L − длина реактора, Pe − число Пекле; U − скорость переноса биомассы, x и z − размерная и безразмерная координаты вдоль оси реактора соответственно.

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

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

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

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

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

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

Устойчивость работы проточных химических реакторов на основе диффузионной модели рассматривалась в ряде работ (напр. [11–14]). Этому фактору уделяется серьезное внимание в литературе. Из предыдущего анализа следует, что при фиксированных значениях параметров 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)(G − b). Использовался метод пристрелки [15]. Для рассматриваемой задачи “стрельбу” целесообразно вести от точки z = 1 до точки z = 0, поскольку при “прямом” направлении интегрирования (от точки z = 0 до точки z = 1) ошибки округления нарастают быстрее, что отмечалось, в частности в работе [16 с. 295]. Для удобства введем новую независимую переменную y = 1 − z и сведем уравнение второго порядка (19) к системе обыкновенных дифференциальных уравнений первого порядка. Пусть T = dG/dy, тогда задача (19), (20) принимает вид

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

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

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

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

Одной из проблем (недостатком) решений задач, полученных сложными, в частности численными алгоритмами, для использования в приложениях является нахождение решения при не приведенных в соответствующих источниках значений параметров. Для инженерной практики обычно привлекательны простые расчетные формулы. Особенно это полезно, когда необходимо использовать функции двух и более переменных. В целях удобной аппроксимации экспериментальных данных в работах [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 выражениями

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

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

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

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

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

ЗАКЛЮЧЕНИЕ

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

ОБОЗНАЧЕНИЯ