Mathematical software of an optimal control system of heat exchange object with dynamic parameters specified in algorithmic form

Full Text

Введение На современных химико-технологических производствах широко распространены аппараты периодического приципа действия, в которых реализуются биохимические и химические процессы, проводимые в условиях натуральной конвекции. Примерами могут служить биореакторы-ферментаторы, применяемые при производстве напитков брожения, дрожжевом производстве, производстве ферментных и медицинских препаратов и т.п. Такие технологические аппараты имеют обычно цилиндрическую или цилиндро-коническую форму и снабжены несколькими рубашками охлаждения (или нагревания, в зависимости от требуемых условий протекания процесса). В промышленности подобные аппараты принято называть танками. Обеспечение свободного конвективного движения жидкой среды в танке (термогравитационной конвекции или конвекции Рэлея) в течение всего периода процесса является зачастую обязательным условием протекания последнего. Так, в пищевой индустрии широкое распространение получили танки значительного объема (50–100 м3) с системой из трех-четырех охлаждающих рубашек, встроенных по вертикальной координате. Отметим, что время протекания (био)технологических процессов в танках обычно значительно и может составлять 10-15 суток [1]. Охлаждение жидкой среды в танке осуществляется с помощью систем автоматического регулирования, установленных для каждой зоны охлаждения, задание регуляторам определяется технологическими условиями процесса. Управляющим воздействием является расход и температура хладагента, подаваемого из холодильных установок. Большая энергоемкость процессов брожения требует применения энергосберегающих алгоритмов управления температурой жидкой среды при гарантированном характере протекания процесса. Другими словами, актуальной является задача достижения минимизации тепловых потерь при гарантированной средней интенсивности теплового потока от жидкой среды в танке к хладагенту в рубашке. Целью данной работы является создание математического обеспечения такого алгоритма. Постановка задачи Будем использовать для математической постановки задачи принципы термодинамики при конечном времени [2], согласно которым между двумя взаимодействующими за заданный (и конечный) временной интервал термодинамическими подсистемами, имеющими разную температуру, возникает тепловой поток, зависящий от термодинамической движущей силы. Диссипацию энергии при теплообмене в этом случае можно определить количественно с помощью производства энтропии , представляющую собой по определению неотрицательную величину и характеризующую непроизводительные потери работоспособной энергии [3]. Практический интерес представляет усредненное за временной интервал протекания процесса производство энтропии , которое для совершенных с точки зрения термодинамики при конечном времени систем должно достигать (или быть максимально приближено) к предельной усредненной оценке , то есть минимально возможному значению. Отметим, что указанная предельная оценка есть величина положительная и обращается в ноль только для идеализированных, обратимых процессов. Разность определяет степень несовершенства системы. Будем рассматривать задачу достижения минимально возможного как оптимизационную с критерием в форме функционала (1) где – независимый аргумент (календарное время); – заданное и конечное время протекания процесса, – усредненная по пространственным координатам фазовая переменная, в рассматриваемом случае – температура жидкой среды в танке. Отметим, что такое усреднение оправдано самой формой критерия; – управление, в рассматриваемом случае – изменяющаяся во времени температура хладагента, находящегося в рубашке. На управление наложено автономное ограничение , то есть , где и – заданные и постоянные величины; – тепловой поток, определяемый характером теплообмена; – термодинамическая движущая сила, обусловливающая интенсивность теплового потока. Отметим, что в соответствии с законом Гюи-Стодолы потери полезной работы, имеющие место в необратимых процессах, прямо пропорциональны приращению энтропии участвующих в процессах теплообмена сред [4]. Поэтому максимальный термодинамический коэффициент полезного действия соответствует минимальному производству энтропии, то есть наименьшему значению задачи (1). Конкретизируем выражения для теплового потока и термодинамической движущей силы. Будем рассматривать ньютоновский теплообмен, характерный для интервалов температур среды и хладагента в описываемых условиях протекания процесса. (2) где – интегральный коэффициент теплопередачи через разделяющую рубашку и танк твердую поверхность. (3) В критериальную постановку (1) необходимо добавить связи и ограничения, обусловливающие характер протекания процесса и требования к его средней интенсивности. Во-первых, введем в рассмотрение требование заданной усредненной во времени интенсивности теплового потока в форме интегрального ограничения (4) где – заданная, постоянная и положительная величина Во-вторых, в постановку задачи необходимо ввести определяющие связи, обусловливающие динамику процесса теплообмена, то есть выражения, описывающие математическую модель процесса. Отметим, что динамика рассматриваемых объектов корректно может быть описана интегро-дифференциальными уравнениями в частных производных, аргументами которых, помимо времени, являются пространственные координаты. Для цилиндрических и цилиндро-конических аппаратов принято применять цилиндрическую систему координат, аргументами которой в общем случае являются: высота слоя среды , радиальная координата и азимутальная координата В инженерных приложениях принято пренебрегать азимутальной координатой и рассматривать три независимых аргумента: время, вертикальную и радиальную составляющую пространственных координат. Кроме того, для процессов, протекающих в условиях естественной термогравитационной конвекции, для плотности среды применяют приближение Буссинеска [5], согласно которому плотность среды считают функцией ее температуры, то есть (5) где: – плотность среды в текущем сечении относительно координаты , – плотность среды при так называемой отсчетной температуре ; в качестве которой в нашем случае может быть выбрана температура среды в нижнем гоизонтальном сечении аппарата, – коэффициент объемного теплового расширения. По Буссинеску, при малой величине и малой вариации параметров среды (кинематической вязкости , температуропроводности и др.) плотность и сами эти параметры считают постоянными во всех членах уравнений, кроме одного: вариация плотности сохраняется там, где она умножается на ускорение силы тяжести [6]. Использование упрощающего приближения Буссинеска позволяет выписать математическую модель объекта в форме системы уравнений Навье-Стокса (6) (7) (8) В выражениях (6) – (8) обозначено: вектор скорости конвективного потока, конвективное ускорение, и вариации температуры и давления соответственно. Символами и = обозначены операторы Гамильтона и Лапласа соответственно. Таким образом, постановка рассматриваемой оптимизационной задачи должна включать в себя критерий (1) и ограничения (4), (6) – (8). Сделанная математическая постановка привела к задаче на принцип максимума Л.С. Понтрягина. Особенностью рассматриваемой задачи является наличие определяющих связей в форме дифференциальных уравнений в частных производных, для которых, как известно, принцип максимума строго не обоснован и традиционный подход к решению применен быть не может. Кроме того, введенные в рассмотрение связи представляют собой выражения, аналитическое решение которых получить чрезвычайно затруднительно, если не невозможно. В таких случаях принято считать, что связи заданы в алгоритмической форме и решение следует искать с помощью численных методов. Откажемся от традиционной схемы решения поставленной задачи. Вместо этого будем использовать двухэтапный алгоритм. На первом этапе отбросим автономные ограничения на управляющую переменную и связи (6) – (8). Таким образом, действуя по логике В.Ф. Кротова, перейдем к расширенной задаче [7], включающей в себя критерий и интегральную связь (4). Аналитическое решение расширенной задачи, то есть оптимальное значение фазовой переменной , будет являться оценочным для исходной постановки. На втором этапе будем осуществлять итерационный поиск наиболее близкого к оценочному физически реализуемого решения. Для невырожденного случая расширенной задачи условие оптимальности примет вид (9) где: – функция Лагранжа, – множитель Лагранжа (постоянный коэффициент). Опуская в обозначениях символы усреднения по пространственным координатам, выпишем для расширенной задачи функцию Лагранжа (10) Необходимые условия оптимальности для (10) в форме требования стационарности по примут вид (11) Отметим, что в поставленной задаче температуры представлены в шкале Кельвина, кроме того, нами принято, что > [0;T]. Учитывая эти обстоятельства, после несложных математических выкладок приходим к соотношению для условий оптимальности (12) Отметим, что в нашем случае Оценочное значение может быть рассчитано из соотношения (13) где константу определим из выражения (14) Второй этап решения. На этом этапе будем считать, что ограничения-связи, определяющие динамику объекта, заданы в алгоритмической форме, то есть для них нет аналитических выражений. Вместо будем получать численное решение модельных уравнений для пошагового расчета допустимых управлений, соответствующих каждому моменту времени и соответствующих им значений критерия оптимальности . Если IД окажется близким к I*, то полученное допустимое решение (допустимое управление) является оптимальным и задача решена. Значительная разница между IД и I* свидетельствует об эффективности расширения. Шаги алгоритма второго этапа решения. 1) Разбиваем время процесса на интервалов , исходя из априорных сведений о динамике процесса. 2) Для первого интервала задаем стартовое значение и рассчитываем по (12) значение управления. Параметр может быть оценен, исходя из условия (15) где: = /, = – безразмерные показатели. Таким образом, стартовое значение может быть выбрано, исходя из желамого (но не обязательно достижимого) соотношения между и . 1) Численно решаем с внутренним шагом систему уравнений (6) – (8) и осуществляем операцию усреднения полученных значений по координатам. При этом для каждого внутреннего шага перед процедурой численного решения пересчитываем и используем его для последующего внутреннего шага. Таким образом, для первого заданного интервала получаем набор табулированных значений и в соответствии с внутренним шагом моделирования. Отметим, что поскольку имеет место итеративный счет, целесообразно применить обычные инженерные допущения, упрощающие вид зависимостей (6) (8). 2) Рассчитываем, используя накопленные в ходе решения модельных уравнений значения и и константу (14). Затем с помощью выражения (13) получаем расчетное значение средней интенсиваности потока . 5) Сравниваем заданное и полученное расчетное значение для средней интенсивности потока: , где – наперед заданное малое положительное число. Если условие не выполнено, то следует изменить значение и повторить шаги 2 4. 6) Если условие шага 5 выполнено или достигнуто другое условие остановки, то полученное итоговое значение следует принять в качестве наиболее близкого к оптимальному оценочному значению и использовать для расчета оптимального допустимого значения управления на заданном временном интервале. На этом этапе следует рассчитать значение критерия оптимальности (одним из методов численного интегрирования) для определения его степени удаленности от критерия оптимальности расширенной задачи, то есть для оценки эффективности расширения. При переходе к следующему временному интервалу полученное на предыдущем шаге оптимальное значение управления и параметр следует использовать в качестве стартовых значений. При синтезе системы организационного управления с использованием предлагаемого алгоритма полученное расчетное значение управления для текущего временного интервала следует использовать в качестве задания регулятору. По истечении реального интервала времени следует скорректировать значение фазовой переменной в соответствии с реальными величинами, полученными от датчика температуры среды. Это скорректированное значение следует использовать в качестве стартового для следующего интервала где описанные выше шаги алгоритма должны быть повторены вплоть до достижения времени окончания процесса . Сделанный авторами вычислительный эксперимент, основанный на реальных данных процесса брожения пива в цилиндро-коническом танке (первый этап брожения), [8] подтвердил работоспособность предлагаемого алгоритма. Заключение Предлагаемый в работе подход позволяет синтезировать систему управления, близкую к термодинамически совершенной. Такая система обеспечивает минимизацию диссипации энергии в среднем за время процесса при гарантированной и заданной интенсивности теплового потока. Практическое применение таких систем особенно актуально для объектов с распределенными параметрами, отличающихся большой энергоемкостью.

About the authors

A. Y. Artyushkin

Moscow State University of Food Production

Email: atp@mgupp.ru
+79037702868

P. V. Voronina

Moscow State University of Food Production

Email: atp@mgupp.ru
Ph.D.; +79037702868

A. V. Tatarinov

Moscow State University of Food Production

Email: atp@mgupp.ru
Ph.D.; +79037702868

References

Supplementary files