Математическое моделирование химико-технологических процессов

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования «НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ ТОМСКИЙ ПОЛИТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ» Н.В. Ушева, О.Е. Мойзес, О.Е. Митянина, Е.А. Кузьменко МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ХИМИКО-ТЕХНОЛОГИЧЕСКИХ ПРОЦЕССОВ Рекомендовано в качестве учебного пособия Редакционно-издательским советом Томского политехнического университета Издательство Томского политехнического университета 2014 УДК 66.011:519.876(075.8) ББК 35.11:22.1я73 М34 Ушева Н.В. М34 Математическое моделирование химико-технологических процессов: учебное пособие / Н.В. Ушева, О.Е. Мойзес, О.Е. Митянина, Е.А. Кузьменко; Томский политехнический университет. − Томск: Изд-во Томского политехнического университета, 2014. – 135 с. В пособии рассмотрена методология построения математических моделей химико-технологических процессов; приведены математические модели структуры потоков, кинетики химических реакций, гомогенных химических реакторов, тепловых и массообменных процессов. Рассмотрены подходы построения математических моделей экспериментально-статистическими методами, методы корреляционного и регрессионного анализа, методы планирования эксперимента и методы оптимизации. Предназначено для студентов, обучающихся по направлениям «Энергои ресурсосберегающие процессы в химической технологии, нефтехимии и биотехнологии», «Химическая технология», «Биотехнология». УДК 66.011:519.876(075.8) ББК 35.11:22.1я73 Рецензенты Доктор технических наук, профессор ТУСУРа С.В. Смирнов Кандидат технических наук заведующая лабораторией Института химии нефти СО РАН Н.В. Юдина © ФГБОУ ВПО НИ ТПУ, 2014 © Ушева Н.В., Мойзес О.Е., Митянина О.Е., Кузьменко Е.А., 2014 © Оформление. Издательство Томского политехнического университета, 2014 ОГЛАВЛЕНИЕ ВВЕДЕНИЕ ............................................................................................................. 6 1. ОБЩИЕ ПРИНЦИПЫ МОДЕЛИРОВАНИЯ .................................................. 8 1.1. Классификация моделей ............................................................................. 8 1.2. Методология построения математических моделей химико-технологических процессов ....................................................... 10 2. ДЕТЕРМИНИРОВАННЫЕ МАТЕМАТИЧЕСКИЕ МОДЕЛИ ХИМИКО-ТЕХНОЛОГИЧЕСКИХ ПРОЦЕССОВ ........................................... 16 2.1. Математическое описание гидродинамической структуры потоков ................................................. 16 2.1.1. Модель идеального смешения ....................................................... 17 2.1.2. Модель идеального вытеснения .................................................... 18 2.1.3. Диффузионные гидродинамические модели ................................ 20 2.1.4. Ячеечные гидродинамические модели ......................................... 21 2.1.5. Определение условий перемешивания в проточных аппаратах ................................................................... 22 2.2. Моделирование тепловых процессов в химической технологии ......................................................................... 24 2.2.1. Основные закономерности теплообмена ...................................... 24 2.2.2. Математические модели теплообменных аппаратов .................. 26 2.2.3. Пример моделирования теплообменных процессов ................... 29 2.3. Математическое моделирование массообменных процессов .............. 31 2.3.1. Математическое описание равновесия в системе «жидкость-пар» и «жидкость-жидкость» .................... 31 2.3.2. Моделирование процесса массопередачи .................................... 35 2.3.3. Моделирование процесса сепарации ............................................ 37 2.3.4. Моделирование процесса ректификации ..................................... 40 2.3.5. Моделирование процесса абсорбции ............................................ 44 2.3.6. Моделирование процесса адсорбции ............................................ 45 2.4. Математическое моделирование кинетики химических реакций ................................................................................. 47 2.4.1. Основные понятия химической кинетики .................................... 47 2.4.2. Моделирование кинетики гомогенных химических реакций .................................................. 50 2.4.3. Моделирование кинетики гетерогенных химических реакций ............................................... 52 2.5. Моделирование гомогенных химических реакторов ............................ 63 2.5.1. Классификация реакторов .............................................................. 64 3 2.5.2. Математическая модель реактора идеального смешения ........... 64 2.5.3. Математическая модель реактора идеального вытеснения ................................................................... 67 2.5.4. Исследование химического процесса, протекающего в гомогенном реакторе идеального смешения ..................................................................... 69 2.5.5. Исследование химического процесса, протекающего в реакторе идеального вытеснения в стационарном режиме .................................................................. 71 3. ЭКСПЕРИМЕНТАЛЬНО-СТАТИСТИЧЕСКИЕ МЕТОДЫ ПОСТРОЕНИЯ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ ........................................ 75 3.1. Основные понятия и определения ........................................................... 75 3.2. Статистические модели объектов на основе пассивного эксперимента ....................................................... 77 3.2.1. Методы корреляционного и регрессионного анализа ................. 78 3.2.1.1. Линейная регрессионная модель с одной независимой переменной .................................... 82 3.2.1.2. Статистический анализ результатов ................................ 83 3.2.1.3. Параболическая регрессионная модель ........................... 86 3.3. Статистические модели на основе активного эксперимента (методы планирования экстремальных экспериментов) ....................... 89 3.3.1. Планы первого порядка .................................................................. 89 3.3.1.1. Полный факторный эксперимент ..................................... 89 3.3.1.2. Статистический анализ уравнения регрессии ................. 95 3.3.1.3. Дробный факторный эксперимент ................................... 96 3.2.3. Статистические модели оптимальной области объекта исследования ............................... 99 3.4. Симплексный метод планирования и оптимизации ............................ 107 4. МЕТОДЫ ОПТИМИЗАЦИИ В ХИМИЧЕСКОЙ ТЕХНОЛОГИИ ........... 112 4.1. Основные понятия и определения ......................................................... 112 4.2. Систематизация методов оптимизации ................................................ 114 4.3. Статистические методы оптимизации .................................................. 115 4.3.1. Метод крутого восхождения по поверхности отклика (Бокса-Уилсона) ............................................................................ 115 4.4. Аналитические методы .......................................................................... 119 4.4.1. Оптимизация реактора идеального смешения ........................... 122 4.4.2. Задача поиска оптимальной температуры обратимой химической реакции .................................................. 124 4.5. Численные методы решения оптимизационных задач без ограничений ............................................ 126 4.5.1. Одномерная оптимизация ............................................................ 126 4.5.1.1. Метод дихотомии ............................................................. 126 4 4.5.1.2. Метод золотого сечения .................................................. 127 4.5.1.2. Метод сканирования ........................................................ 129 4.5.2. Многомерный поиск оптимума ................................................... 130 СПИСОК ЛИТЕРАТУРЫ .................................................................................. 133 5 ВВЕДЕНИЕ Математическое моделирование – метод исследования процессов или явлений на математических моделях с применением ЭВМ. Современный уровень развития вычислительной техники расширяет возможности использования метода математического моделирования при исследовании кинетики гомогенных и гетерогенных химических реакций, лежащих в основе промышленных процессов; выборе типа химического реактора, теплообменного и массообменного оборудования; получении оперативных прогнозов и решении задач оптимизации технологических режимов ведения промышленных процессов действующих производств в условиях меняющихся состава сырья и производительности, а также при проектирования технологических схем новых и модернизируемых производств химической промышленности. Процессы, связанные с химической технологией, очень сложны. Это прежде всего химические превращения в аппаратах различных конструкций, обусловленных особенностями протекания химических реакций, многокомпонентностью и многостадийностью многих из них, необходимостью проведения катализа. Не менее сложны и массообменные процессы, в частности процессы ректификации многокомпонентных смесей, широко используемые при подготовке сырья для химических превращений и разделении продуктов реакций либо отделения непрореагировавших компонентов сырья от продуктового потока. В настоящее время широко используются совмещенные реакционно-ректификационные процессы как более энерго- и ресурсосберегающие и эргономичные. Теплообменные процессы являются неотъемлемой частью любого химического производства. Их эффективность зависит от конструкций аппаратов, свойств теплоносителей и ряда технологических параметров. Поэтому важным этапом математического моделирования является создание математической модели, которая бы адекватно описывала рассматриваемый процесс. Обычно создаются математические модели отдельных аппаратов, базирующиеся на моделях процессов, протекающих в этих аппаратах, а затем моделируются технологические схемы, связывающие эти аппараты в единый технологический процесс. В зависимости от сложности самого процесса и возможностей получения экспериментальной информации о его прохождении, при разработке математических моделей используется либо детерминированный подход, в основе которого лежат фундаментальные законы, либо 6 эмпирический, в основе которого лежит статистическая обработка экспериментальной информации. Поскольку математические модели могут быть представлены линейными, нелинейными, дифференциальными уравнениями, уравнениями в частных производных и их системами, в зависимости от сложности моделируемых явлений, необходимо знать и уметь применять численные методы для их решения. Чтобы решение задач оптимизации было реализуемо, нужно правильно определить критерии оптимальности, представить функцию цели, задать ограничения на оптимизирующие параметры и грамотно выбрать метод оптимизации. И наконец, чтобы воспользоваться вычислительной техникой и решить уникальную задачу, связанную с моделированием конкретного химико-технологического процесса, необходимо знать какой-либо из современных языков программирования и уметь работать в соответствующей среде, создавая удобный для пользователя интерфейс. Конечно, для решения задач выбора наиболее подходящего численного метода могут быть привлечены математики, для создания программы с удобным для пользователя интерфейсом – профессиональные программисты, но саму математическую модель должны создавать специалисты предметной области, т. е. специалисты, компетентные в области химической технологии и промышленной реализации химических и нефтехимических производств. 7 1. ОБЩИЕ ПРИНЦИПЫ МОДЕЛИРОВАНИЯ 1.1. Классификация моделей В настоящее время моделирование широко используется в различных областях науки и техники. Широкое применение моделей объясняется тем, что модель дает возможность установить в явлении, объекте или процессе основные закономерности, которые им присущи, и пренебречь второстепенными, вспомогательными признаками [1, 2]. На рис. 1.1 приведена общая классификация моделей. Рис. 1.1. Общая классификация моделей В зависимости от характера и сложности явлений могут использоваться различные методы моделирования:  геометрический (на основе геометрического подобия величин);  физический (характеризуется одинаковой физической природой модели и исследуемого объекта);  математический (характеризуется различной физической природой и одинаковым математическим описанием модели и исследуемого объекта). Процессы химической технологии – это сложные физико-химические системы, имеющие двойственную детерминировано-стохастическую природу, переменные в пространстве и во времени. Особенности данных процессов состоят в следующем [3]:  в участии многокомпонентных и многофазных материальных потоков; 8  наличии процессов переноса импульса, энергии, массы на границе раздела фаз;  на процесс в значительной степени влияют геометрические характеристики аппарата;  наложении стохастических особенностей гидродинамической обстановки в аппарате на процессы массо-, теплопереноса и химического превращения. Это объясняется случайным взаимодействием составляющих компонентов фаз (соударением частиц, коалесценцией) или случайным характером геометрии граничных условий в аппарате. Подобного рода системы характеризуются чрезвычайно сложным взаимодействием составляющих их фаз и компонентов, вследствие чего изучение их с позиции классических детерминированных законов переноса и сохранения становится невозможным. Ключ к решению данной задачи дает применение метода математического моделирования, базирующегося на основе стратегии системного анализа, сущность которого заключается в представлении процесса как сложной взаимодействующей иерархической системы с последующим качественным анализом ее структуры, разработкой математического описания и оценкой неизвестных параметров [3]. Математическим моделированием называют изучение свойств объекта на математической модели, целью которого является определение оптимальных условий протекания процесса, управление им на основе математической модели и перенос результатов на объект [1–4]. Математическая модель химико-технологического процесса (ХТП) – совокупность математических структур: формул, уравнений, неравенств и т. д., адекватно описывающая исследуемые свойства объекта. Реализованная на компьютере математическая модель называется компьютерной математической моделью, а проведение целенаправленных расчетов с помощью компьютерной модели называется вычислительным экспериментом. Математическое моделирование включает в себя три взаимосвязанных этапа:  составление математического описания изучаемого объекта. Применительно к химической технологии математическая модель – совокупность математических зависимостей, отражающих в явной форме сущность химико-технологического процесса и связывающих его физические, режимные, физико-химические и конструктивные параметры;  выбор метода решения системы уравнений математического описания и реализация его в форме моделирующей программы;  установление соответствия (адекватности модели объекту). 9 В модели должны быть учтены все наиболее существенные факторы, влияющие на процесс, и вместе с тем она не должна быть загромождена множеством мелких, второстепенных факторов, учет которых только усложнит математический анализ. В зависимости от конкретной реализации процесса и его аппаратурного оформления, все многообразие химико-технологических процессов можно разделить на четыре класса:  процессы, переменные во времени (нестационарные);  процессы, не меняющиеся во времени (стационарные);  процессы, в ходе которых их параметры не изменяются в пространстве;  процессы с учетом пространственного изменения параметров. Так как математические модели являются отражением соответствующих объектов, то они классифицируются аналогичным образом [1–4]. 1.2. Методология построения математических моделей химико-технологических процессов В зависимости от подхода к формированию математического описания и природы процессов, протекающих в моделируемых объектах, различают два класса моделей: стохастические и детерминированные [1, 3, 4]. Стохастические (эмпирические, статистические) модели – отражают вероятностный характер явлений, когда рассчитывается не истинное значение параметров процесса, а вероятность их расчета в определенном интервале значений. Данные модели не несут информации о физикохимической сущности решаемой задачи, но их простота позволяет их эффективно использовать при моделировании химико-технологических процессов (ХТП). Стохастическая модель описывает процесс, в котором значение выходной величины не находится в однозначном соответствии с входной величиной. Пример: формула Войнова для расчета молекулярной массы узких нефтяных фракций по их средней температуре кипения М = 52,63 + 0,246Т + 0,01Т2. Детерминированные (причинные, структурные, знаковые) модели отражают детерминированную (причинную) сущность взаимосвязи исследуемых явлений, когда можно теоретически обосновать изменение поведения системы; объясняют сущность взаимосвязи явлений, протекающих в моделируемой системе и описываемых уравнениями статики и динамики химических, физико-химических, тепловых, гидродинамических процессов химической технологии. 10 Детерминированная модель описывает процесс, в котором значение выходной величины однозначно определяется значением входной величины В качестве примера детерминированной модели можно привести уравнение Аррениуса, описывающее влияние температуры T на величину константы скорости химической реакции k, справедливое для любых реакций: k  k0  e E RT , где E – энергия активации; R – универсальная газовая постоянная; k0 – предэкспоненциальный множитель. Математическая модель является системой уравнений математического описания, отражающей сущность протекающих в объекте явлений, для которой определен алгоритм решения, реализованный в форме моделирующей программы. Согласно этому определению, математическая модель должна рассматриваться в совокупности трех ее аспектов: смыслового, аналитического, вычислительного [3]. Основные аспекты математической модели представлены на рис. 1.2. Рис. 1.2. Составные части математической модели 11 Физическое описание природы моделируемого объекта Построение любой модели начинается с физического описания объекта моделирования. При этом выделяют «элементарные» процессы, протекающие в объекте моделирования, которые подлежат отражению в модели [3, 4]. Под элементарным процессом в данном случае понимается физико-химический процесс, относящийся к определенному классу явлений. Обычно при моделировании объектов химической технологии учитывают следующие «элементарные» процессы:  движение потоков фаз;  массообмен между фазами;  теплопередачу;  изменение агрегатного состояния (испарение, конденсация, растворение и т. д.);  химические превращения. На практике часто делают допущения относительно характера связей между «элементарными» процессами с целью избежать излишнего усложнения в описании. Составление математического описания объекта При составлении математического описания процесса, как правило, используют блочный принцип (системный подход), согласно которому составлению математического описания предшествует анализ отдельных «элементарных» процессов, протекающих в объекте моделирования. При этом эксперименты по изучению каждого такого процесса проводят в условиях, максимально приближающихся к условиям эксплуатации объекта моделирования. Далее составляется математическое описание каждого из этих процессов. Заключительным этапом является объединение всех исследованных элементарных процессов (блоков) в одну систему уравнений математического описания объекта моделирования. Достоинством данного принципа является то, что его можно использовать на стадии проектирования объекта, когда окончательный вариант аппаратурного оформления еще неизвестен. К методам составления математического описания ХТП относят аналитический, экспериментальный и экспериментально-аналитический. Характеристика данных методов приведена на рис. 1.3. В составе математического описания на основе физической природы объекта можно выделить несколько групп уравнений. Формирование математической модели в виде совокупности подсистем (блоков) приведено на рис. 1.4. Применение метода математического моделирования связано с большим объемом расчетов по модели, реализуемым на компьютерах, что, в свою очередь, вызвано большим числом уравнений, входящих в математическое описание задачи [6]. 12 13 Рис. 1.3. Методы формирования математической модели Рис. 1.4. Математическое описание химико-технологического процесса Алгоритм моделирования ХТП включает несколько взаимосвязанных этапов: 1. Формирование исходных данных моделирования. На основе литературных данных собирается, перерабатывается необходимая информация о моделируемом процессе: формируется список уравнений математического описания, банк данных, необходимый для расчета уравнений математического описания, анализируется возможное аппаратурное оформление процесса. На этом этапе формируются детерминированные математические модели. 2. Формирование математической модели ХТП. Осуществляется разработка модели с последовательным переходом от низшего уровня иерархии моделируемого объекта к высшему – от единичного элементарного акта к реализации процесса в конкретном аппаратурном оформлении. На данном этапе также формируются начальные и граничные условия реализации процесса. 3. Корректное упрощение математической модели. Осуществляется поиск путей упрощения математической модели с целью сокращения времени расчета без значительного искажения результатов последующего моделирования, например: допущение о квазигомогенности среды (в этом случае нет необходимости учитывать процессы, протекающие на границе раздела фаз), допущение об изотермичности режима работы реактора (при этом не требуется знание зависимости констант скорости химической реакции от температуры) и т. д. 4. Выбор алгоритма решения математической модели. Решение модели требует применения численных методов расчета, при этом одна 14 и та же задача может быть решена различными методами, отличающимися друг от друга сложностью программирования алгоритма и быстродействием. Например, для решения алгебраического уравнения можно использовать методы половинного деления, сканирования, касательных, секущих и т. д. 5. Разработка программы расчета, отладка, получение пробных решений. Один из наиболее трудоемких этапов компьютерного моделирования. 6. Оценка адекватности разработанной математической модели. Выполняется сопоставление результатов расчета по модели с практическими данными, полученными в ходе экспериментов на реальном объекте. 7. Интерпретация результатов вычислительного эксперимента и выдача практических рекомендаций. Вопросы для самоконтроля 1. 2. 3. 4. 5. 6. Дайте определение понятиям «модель», «моделирование». Приведите примеры. Перечислите виды моделирования, проанализируйте возможности их применения в химической технологии. В чем заключается основная сложность моделирования химико-технологических процессов? Назовите два основных вида математических моделей. Приведите примеры. В чем отличие стохастических моделей от детерминированных? Перечислите основные этапы математического моделирования. К каким этапам моделирования необходимо вернуться, если расчет на модели показал неадекватный результат? 15 2. ДЕТЕРМИНИРОВАННЫЕ МАТЕМАТИЧЕСКИЕ МОДЕЛИ ХИМИКО-ТЕХНОЛОГИЧЕСКИХ ПРОЦЕССОВ 2.1. Математическое описание гидродинамической структуры потоков Химико-технологические процессы обычно протекают в движущихся потоках, гидродинамические закономерности перемещения которых оказывают существенное влияние на эффективность химических производств, поэтому при составлении математических моделей ХТП важное значение приобретает описание движения потоков веществ. Движущиеся потоки могут быть как одно- так и многофазными и обычно имеют сложную структуру. В гидродинамике существует целый ряд уравнений, с помощью которых можно описать движение среды (например, уравнение Навье-Стокса, уравнение неразрывности потока и др.) [7]. Применение классических законов гидродинамики к химикотехнологическим процессам оказывается затруднительным из-за сложности уравнений гидродинамики реальных потоков, поэтому на практике при составлении математических моделей гидродинамических потоков обычно используют более простые приближенные представления об их внутренней структуре. Структура движущейся технологической среды характеризуется степенью перемешивания частиц потока, которая определяет поле концентраций и градиенты температур. Это и послужило предпосылкой для установления по признаку перемешивания некоторых типовых моделей движущихся потоков. Принципиально можно построить гидродинамические модели потоков различной сложности, наиболее отвечающие применяемым конструкциям технологического оборудования. Обычно при составлении математических моделей ХТП используют приближенные (модельные) представления о структуре движущихся потоков отдельных фаз. Это следующие гидродинамические модели:  идеального смешения;  идеального вытеснения;  диффузионные (одно- и двухпараметрические) модели;  ячеечные модели;  комбинированные модели. 16 При построении модели структуры потока должны учитываться следующие требования:  модель должна отражать физическую сущность реального потока и при этом должна иметь достаточно простое математическое описание;  должна давать возможность определять ее параметры расчетным или экспериментальным способом;  должна быть удобной для использования при расчетах конкретного ХТП. 2.1.1. Модель идеального смешения По данной модели поток представляется в виде непрерывной среды, которая поступает в аппарат и мгновенно распределяется по всему объему аппарата вследствие полного (идеального) перемешивания частиц потока, при этом концентрация и температура остаются постоянными во всех точках объема данного аппарата и на выходе из него [4]. Свх v вх С Свых v вых Cвх  С0 ,  вх   0 при постоянном объеме ( V  const ). Уравнение материального баланса потоков на входе и выходе из аппарата: I вх   Свх , I вых   Свых , где I – поток вещества [моль/с], v – объемный расход потока, м3/с; Свх, Свых, С – концентрация вещества в потоке на входе в аппарат, выходе из него и в любой точке объема аппарата соответственно, моль/м3; V – объем, м3. В установившемся режиме I вх  I вых . Если на входе в аппарат произошло изменение концентрации (возмущение), то I вх  I вых и в аппарате произойдет накопление вещества. Предположим, что рассматриваемое 17 изменение в аппарате произошло за очень маленький промежуток времени t  dt , за который в аппарате произойдет накопление массы: ∆М →dM. Разделив обе части уравнения на объем аппарата (V), получим М  d   V     C  C , M  C;  0  (2.1) dt V V dC     C0  C  . dt V Уравнение (2.1) описывает изменение концентрации в аппарате идеального смешения. Учитывая, что время контакта равно  V / v , получим модель идеального смешения в следующем виде: dC 1    C0  C  . (2.2) dt  Начальные условия: при t  0 C  0   C0 . Гидродинамическая модель идеального смешения является моделью с сосредоточенными параметрами, т. к. переменная С изменяется только во времени. Модель идеального смешения (МИС) обычно используют при описании аппаратов, в которых обеспечивается интенсивное перемешивание сред. Это аппараты небольших размеров с соизмеримыми высотой и диаметром. На практике – это аппараты с мешалками барботажного типа либо аппараты с очень высокой скоростью циркуляции потока. 2.1.2. Модель идеального вытеснения Модель идеального вытеснения – это модель с идеализированной структурой потока, в котором принимается поршневое течение без перемешивания вдоль потока при равномерном распределении концентраций вещества в направлении, перпендикулярном движению потока [4]: Свх Cвых uвх uвых Здесь uвх, uвых – линейная скорость потока на входе и выходе из аппарата, м/с. 18 Выделим некоторую элементарную ячейку: Ij–1 Cj–1 Ij Cj Поток на входе в j–1-е сечение равен I вх  u  S  C j 1, где S – площадь поперечного сечения аппарата, м2. Поток на выходе из j-го сечения I вых  u  S  C j . В установившемся режиме I вх  I вых . При изменении концентрации на входе в аппарат I вх  I вых . В системе за некоторый промежуток времени t происходит накопление вещества (dM). Считаем, что при t  0 V  dV , t  dt . Тогда dM   I вх  t   I вых  t   dt   I j 1  t   I j  t    dt; dM    I j  t   I j 1  t   dt; dM   uSC j (t )  uSC j 1 (t )   dt ; dM  uSdC. dt Разделив обе части уравнения (2.3) на dV, получим C t  uS C V (2.3) . Таким образом, уравнение гидродинамической модели идеального вытеснения (МИВ) будет иметь вид C t  u  C l , (2.4) где u – линейная скорость потока, м/с; l – длина аппарата, м; t – время, с. Начальные условия: при t = 0 С(0, l) = C0; граничные условия: при l = 0 C(t, 0) = С0. 19 Это модель с распределенными параметрами, т. к. концентрация вещества в потоке изменяется во времени и по длине аппарата. На практике встречается достаточно большое многообразие аппаратов, которые можно описать моделью идеального вытеснения: химические реакторы (трубчатого, полочного типа и др.), массообменные и теплообменные аппараты. Известно, что если отношение длины (высоты) аппарата к его диаметру ≥20, то в таком аппарате возможна реализация режима идеального вытеснения. 2.1.3. Диффузионные гидродинамические модели Диффузионные модели применяют при описании реальных потоков в аппаратах, в которых происходит продольное и радиальное перемешивание веществ. Природа возникновения продольного и радиального перемешивания весьма сложна. Предположим, что перемешивание в потоке идеального вытеснения (конвективный поток) возникает за счет молекулярной диффузии. Тогда в соответствии с блочным принципом построения математических моделей уравнение однопараметрической диффузионной гидродинамической модели примет следующий вид [4, 5]: C C  2C  u  DL , (2.5)  t  l  l2 где DL – коэффициент диффузии в продольном направлении, м2/с; Начальные и граничные условия: при t = 0 C(0,l) = C0; при l = 0 dC/dt = 0. Однопараметрическая модель учитывает диффузию только в продольном направлении. Если в потоке одновременно происходит продольное и радиальное перемешивание веществ, то математическое описание гидродинамики потока можно представить двухпараметрической диффузионной моделью:   2C 1  C   C  C  2C  U  DL  DR   (2.6) , 2  t  l  l2 r r  r где DR – коэффициент диффузии в продольном направлении, м2/с; r – текущий радиус, м; R – радиус аппарата, м. Для оценки степени влияния диффузии можно использовать диффузионный критерий Пекле [7] U l PeD  . DL При PeD >200 движущийся поток можно считать потоком идеального вытеснения. 20 Диффузионные модели достаточно точно отражают структуру потоков во многих реальных аппаратах: пленочных, распылительных, барботажных колоннах, экстракторах и др. 2.1.4. Ячеечные гидродинамические модели Физическая сущность ячеечной модели заключается в том, что материальный поток может быть представлен несколькими последовательно соединенными ячейками, при этом допускается, что в каждой ячейке поток имеет структуру идеального смешения, а между ячейками перемешивание отсутствует. C1 Свх, v вх C1 C2 C n–1 Ci .... C1 Cn–1 Ci–1 C2 ... Ci Cn, Пусть Vi – объем i-й ячейки, м3. Примем V1 = V2 = V3 = … = Vn. Поскольку в каждой ячейке реализуется режим идеального смешения, то для любой ячейки справедливо уравнение МИС:  dC1 1    C0  C1  ;  dt 1  .  .  dCi 1    Ci 1  Ci  ; (2.7)   dt i  .  .  dC 1  n    Cn1  Cn  .  dt  n Время пребывания вещества в каждой ячейке τ = Vi/v, общее время пребывания τ = V/v, тогда объем всех ячеек V = nVi, где n – число ячеек. 21 Уравнение ячеечной модели для i-й ячейки примет вид dCi n    Ci 1  Ci  , (2.8) dt  при t = 0 C(0) = C0. При n = 1 получим модель идеального смешения, при n →∞ – модель идеального вытеснения. При использовании ячеечной модели очень важно правильно выбрать число ячеек, которое можно рассчитать по формуле Pe ul n D  . 2 2 DL Ячеечные модели достаточно точно отражают свойства потоков в различных абсорбционных и экстракционных колоннах, в теплообменных аппаратах некоторых конструкций, в каскаде химических реакторов с мешалками, в аппаратах с псевдоожиженным слоем. 2.1.5. Определение условий перемешивания в проточных аппаратах Для того чтобы установить характер перемешивания потока в аппарате, необходимо на входе в поток ввести какое-либо вещество (индикатор, трассёр) и изучить изменение концентрации этого вещества в выходном потоке, т. е. найти отклик системы на входное возмущение. Индикатором является вещество (азот, аргон, гелий), которое вводится в небольшом количестве и отличается по свойствам от других компонентов потока. N2 СО+Н2 CH3OH Для измерения концентрации на выходе из аппарата можно использовать различные физико-химические методы анализа: хроматографические, спектральные и др. 22 Существует несколько стандартных способов ввода индикатора в поток:  импульсный;  ступенчатый. При импульсном вводе индикатор вводится в основной поток за короткий промежуток времени и в небольшом объёме, затем снимается изменение концентрации индикатора во времени на выходе из аппарата, т. е. получают кривую отклика системы на импульсное возмущение (С – кривая). Ступенчатый ввод индикатора предполагает замену части основного потока индикатором, при этом на выходе получаем кривую отклика, которая называется F-кривая. В табл. 2.1 приведены кривые отклика для различных гидродинамических моделей. 1. 2. 3. 4. 5. 6. 7. Вопросы для самоконтроля Назовите типовые математические модели структуры потоков в аппаратах. Что такое кривая отклика? Перечислите методы определения гидродинамической структуры потоков. Перечислите модели идеального вытеснения. Перечислите модели идеального смешения? Дать характеристику диффузионной модели? Дать характеристику ячеечной модели? Таблица 2.1 Характер отклика при возмущении: ступенчатом импульсном Модель F(t) С(t) Идеального вытеснения  F( t) С(t) Идеального смешения 23 Окончание табл. 2.1 Характер отклика при возмущении: ступенчатом импульсном Модель С(t) F(t) Диффузионная F(t) С(t) Ячеечная  2.2. Моделирование тепловых процессов в химической технологии 2.2.1. Основные закономерности теплообмена Тепловые процессы в химической технологии имеют как самостоятельное значение при сушке, выпаривании, нагревании, охлаждении и т. д., так и сопровождают химические и массообменные процессы [1, 8]. На рис. 2.1 приведены примеры теплообменных аппаратов. Теплообмен обусловлен стремлением системы к тепловому равновесию. Связь между градиентом температуры и молекулярным потоком теплоты (qT, Вт/м2) определена законом теплопроводности Фурье [7] q    gradT , T где  – коэффициент теплопроводности среды, Вт/(мК). При движении в жидкостях и газах происходит конвективный перенос энергии веществом: qk    I  u , где u – скорость движения среды, м/с;  – плотность вещества, кг/м3; I – энтальпия, Дж/кг. 24 а б Рис. 2.1. Теплообменные аппараты: а – пластинчатый; б – кожухотрубный При конвективном теплообмене плотность теплового потока q определяется суммой молекулярной и конвективной составляющих: q  qT  qk    gradT    I  u . Для упрощения расчета переноса теплоты между поверхностью твердого тела и движущейся сплошной средой используют закон теплоотдачи Ньютона-Рихмана Q    Tc  Tср   F , где  – коэффициент теплоотдачи, Вт/(м2К); F – поверхность теплообмена, м2; Тс – температура стенки, К; Тср – температура среды, К. Коэффициент теплоотдачи зависит от скорости движения жидкости, ее плотности и вязкости, от тепловых свойств жидкости (удельной теплоемкости, теплопроводности), от формы и определяющих размеров стенки и других факторов. Теплоотдача определяется не только тепловыми, но и гидродинамическими условиями. Поэтому конвективный теплообмен описывается дифференциальным уравнением Фурье-Кирхгофа [7]   2T  2T  2T  Т T T T u u u  a  2   , xx y y z z  y 2  z 2  t x где a   / C p   – коэффициент температуропроводности, м2/с; t – время, с; Ср – теплоемкость, Дж/кг К. Количество тепла, передаваемое от нагретого теплоносителя к холодному, определяется основным уравнением теплопередачи [8, 10] Q  К Т  F  T , 25 где КТ – коэффициент теплопередачи, Вт/(м2С); ∆Т– средняя разность температур между теплоносителями. При теплопередаче через стенку толщиной с коэффициент теплопередачи можно рассчитать с помощью уравнения аддитивности термических сопротивлений на пути теплового потока 1 1  1   с  r31  r32  , 2 KТ 1 с где 1 и 2 – коэффициенты теплоотдачи от жидкости к стенке и от стенки к другой жидкости соответственно, Вт/(м2К); с – теплопроводность материала стенки, Вт/(мК); rЗ1 и rЗ2 – термические сопротивления слоев загрязнений с обеих сторон стенки, м2К/Вт. Это уравнение справедливо для передачи тепла через плоскую или цилиндрическую стенку при условии, что RH RB  2 ( Rн и Rв – наружный и внутренний радиусы цилиндра соответственно). 2.2.2. Математические модели теплообменных аппаратов При построении математических моделей теплообменных аппаратов предварительно проводится анализ структуры движения потоков в аппарате. Для каждого потока записывается математическое описание в виде выражения, характеризующего изменения температуры в потоке теплоносителя во времени, обусловленное движением потока и теплопередачей [1, 8]. Предварительно формулируются допущения. Если структура потока теплоносителя соответствует модели идеального смешения, то математическое описание этого потока, с учетом теплопередачи, будет иметь вид dT  vC p Tвх  T   KT F T , V C p (2.9) dt где V – объем потока идеального смешения, м3;  – плотность теплоносителя, кг/м3; Ср – удельная теплоемкость теплоносителя, Дж/(кгК); v – объемная скорость потока, м3/с; F – поверхность теплообмена, м2; КТ – коэффициент теплопередачи, Вт/(м2К); ∆Т– средняя разность температур между теплоносителями; Твх – температура потока на входе, К; t – время, с. Если структура потока соответствует модели идеального вытеснения, то с учетом теплопередачи можно записать: T Т F      C p  K F T , S    Cp (2.10)  t  l L T где S – площадь поперечного сечения потока, м2; L – длина зоны идеального вытеснения, м; l – пространственная координата, изменяющаяся 26 от 0 до L; Т = Т (l, t) – функция распределения температуры потока теплоносителя по пространственной координате во времени. Обычно в уравнениях (2.9) и (2.10) принимают коэффициент теплоотдачи, плотность и теплоемкость теплоносителя постоянными в исследуемом интервале изменения температуры. Предполагается, что объемные скорости потоков остаются постоянными. Рассмотрим математические модели некоторых типов теплообменных аппаратов. Теплообменник типа «смешение-смешение». Примем, что тепло передается от первого потока теплоносителя ко второму. Режим движения потоков – идеальное смешение (рис. 2.1,а). Рис. 2.1,а Схематическое изображение теплообменника типа «смешение-смешение»  Если тепловой емкостью стенки, разделяющей потоки теплоноси- телей, можно пренебречь, то математическая модель аппарата будет состоять из двух уравнений теплового баланса: dT1  V C    1  1  C p1  T1,0  T1   F1  K T  T1  T2  ;  p 1 1 1  dt (2.11)  dT 2 V    C  2   2  C p 2  T2,0  T2   F2  K T  T1  T2  .  2 2 p 2 dt  Если тепловой емкостью стенки, разделяющей потоки теплоносителя, пренебречь нельзя, то необходимо к уравнениям (2.11) добавить уравнение изменения температуры стенки: dT1   V C    1  1  C p1  T1,0  T1   F1  1  T1  T3  ; p 1 1 1  dt   (2.12) dT2      ; V C C T T F T T                p2 2 2 2,0 2 2 2 3 2  2 2 p 2 dt  G  C dT3  F    T  T   F    T  T  , 1 1 1 3 2 2 3 2  3 3 dt 27 где G3 – масса материала стенки, кг; С3 – удельная теплоемкость матеДж риала стенки, ; Т3 – температура стенки, К; 1, 2 – коэффициенты кг  К теплоотдачи, Вт/м2К. Теплообменник типа «смешение – вытеснение» (рис. 2.2) v1 , Т v1, Т 10 Т3 v ,Т 2 20 1 v ,Т 2 2 Рис. 2.2. Схематическое изображение теплообменника типа «смешение – вытеснение»  Тепловой баланс без учета теплоемкости стенки: dT1  1  1  C p1  Т1,0  Т1   F1  K T  T1  T2  ; dt (2.13)  T2  Т 2 F2  2   2  C p 2   K T  T1  T2  . S2  2  C p 2 t l L  С учетом теплоемкости стенки dT V1  1  C p1 1  1  1  C p1  Т1,0  Т1   F  1  T1  T3  ; dt T G3  C3 3  F1  1  T1  T3   F2   2  T3  T2  ; (2.14) t T Т F S2  2  C p 2 2  2  2  C p 2 2    2  T3  T2  . t l L Теплообменник типа «вытеснение – вытеснение» (рис. 2.3). V1  1  C p1 1, Т 10 1, Т 1 Т3 2,Т 20 2, Т 2 Рис. 2.3. Схематическое изображение теплообменника типа «вытеснение-вытеснение» 28  Уравненияе теплового баланса без учета тепловой емкости стенки:  T1 Т F  1  1  C p1 1  1  K T  T1  T2  ; t l L T Т F S 2   2  C p 2 2  v2   2  C p 2 2  2  KT  T1  T2  . t l L S1  1  C р1 (2.15)  С учетом теплоемкости стенки  T1 Т F  1  1  C p1 1  1  1  T1  T3  ; t l L T Т F S2  2  C р 2 2  2  2  C p 2 2  2   2  T3  T2  ; t l L T G3  C3 3  F1  1  T1  T3   F2   2  T3  T2  . t S1  1  C p1 (2.16) 2.2.3. Пример моделирования теплообменных процессов В теплообменнике типа «труба в трубе» охлаждается жидкость. Хладоагент и охлаждающаяся жидкость движутся прямотоком. Необходимо рассчитать температуру теплоносителей на выходе из аппарата и получить температурные профили по длине аппарата. В табл. 2.2 приведены исходные данные для расчета. Таблица 2.2 Теплоноситель горячий холодный 200 35 –4 5,110–4 2,310 900 1000 3 3,3510 4,19103 0,01 0,03 Параметры Температура, °С Объемная скорость, м3/с Плотность, кг/м3 Теплоемкость, Дж/кг°С Диаметр трубы, м В теплообменнике реализуется режим «вытеснение-вытеснение», поэтому математическое описание будет иметь вид T Т F S1  1  C p1 1  1  1  C p1 1  1  KT  T1  T2  ; t l L T Т F S2  2  C p 2 2  2  2  C p 2 2  2  KT  T1  T2  . t l L В стационарном режиме работы теплообменника, когда ∂T1 /∂t = 0; ∂T2 /∂t = 0, уравнения теплового баланса примут следующий вид: 29 KT    d  dT1    T1  T2  ;  dl     C p1 1 1    dT2  KT    d  T  T  ,  dl 2   2  C p 2 1 2 где d – диаметр трубы теплообменника, м. Для удобства вычисления введем обозначения: K   d ; b1  1  1  C p1 b2  (2.17) K   d . 2   2  C p 2 Систему полученных дифференциальных уравнений (2.17) решаем с помощью численного метода Эйлера [12]: T i  T i1  hb T i1  T i1 ; 1 1 1 2  1  i T2  T2i1  hb2 T1i 1  T2i1 ,  где i – номер шага по длине теплообменника; h – шаг интегрирования по длине теплообменного аппарата. На рис. 2.4 приведены результаты расчета процесса теплообмена. С применением данной математической модели можно выполнить исследования влияния температуры, расхода теплоносителя и хладоагента, размеров аппарата на процесс теплообмена.     200 Температура, ◦С 180 160 140 120 Т1 100 80 60 Т2 40 20 0,5 1 1,5 2 Длина теплообменника, м 2,5 Рис. 2.4. Изменение температуры по длине теплообменного аппарата 30 3 3,5 Вопросы для самоконтроля 1. Перечислите основные тепловые процессы в химической технологии. 2. Какие гидродинамические модели структуры потоков применяются при моделировании теплообменных аппаратов? 3. Перечислите параметры математической модели теплообменных аппаратов и их размерности. 4. Каковы принципы составления уравнений тепловых балансов? 5. Перечислите управляющие параметры процесса теплообмена. 6. В чем отличие математической модели трубчатой печи от модели теплообменного аппарата? 7. На основании каких законов разрабатываются математические модели тепловых процессов? 8. Дать характеристику математической модели теплообменного аппарата типа «смешение-смешение». 9. Дать характеристику математической модели теплообменного аппарата типа «вытеснение-вытеснение». 10. Дать характеристику математической модели теплообменного аппарата типа «перемешивание-вытеснение». 2.3. Математическое моделирование массообменных процессов Любой массообменный процесс представляет собой сложную систему, состоящую из ряда подсистем: «равновесие», «массопередача», «гидродинамика», «теплопередача», «подсистема балансов массы и энергии» [12, 13]. При этом математическое описание массообменного процесса создается на основе математических моделей отдельных блоков, каждая из которых характеризуется собственным набором параметров. Рассмотрим подробнее моделирование равновесия в массообменных системах. 2.3.1. Математическое описание равновесия в системе «жидкость-пар» и «жидкость-жидкость» Равновесие в системе «жидкость-пар» При расчете равновесия «жидкость-пар» можно выделить четыре типа задач, в зависимости от того, какие переменные задаются и какие рассчитываются:  расчет состава пара и температуры смеси по известному составу жидкости и давлению;  расчет состава пара и давления по составу жидкости и температуре;  определение состава жидкости по составу пара при известном давлении;  определение состава жидкости по составу пара при известной температуре. 31 Для каждого компонента i в смеси условие термодинамического равновесия в системе задается уравнением [12] fi G  fi L , (2.18) где f – фугитивность; индексы G и L относятся к паровой и жидкой фазе соответственно. Фундаментальной задачей моделирования является установление связей фугитивностей с составами смесей, поскольку при разработке химико-технологических процессов используются именно эти составы. Фугитивность компонента в смеси зависит от состава смеси, температуры и давления. Для связи fiG с температурой, давлением и составом используют коэффициент фугитивности [14] Фi  fi G . P  yi (2.19) Здесь yi – мольная доля компонента в паровой фазе; P – давление в системе, Па. Для смеси идеальных газов Фi = 1. Фугитивность компонента i в жидкой фазе связана с составом фазы и коэффициентом активности соотношением i  i xi  fi L , xi  fi 0 (2.20) в котором xi – мольная доля компонента в жидкой фазе; γi – коэффициент активности компонента i; fi0 – стандартная фугитивность компонента i (фугитивность чистой жидкости i при общем давлении в системе P и xi = 1). Классическая термодинамика не позволяет определить вид зависимости коэффициента активности от состава и температуры. Однако существует термодинамическое соотношение, позволяющее коррелировать и обобщать экспериментальные данные – уравнение ГиббсаДюгема, согласно которому коэффициенты активности связаны между собой следующим соотношением [11]: N   ln  i    0.  xi T , P  x  i 1 i (2.21) С практической точки зрения уравнение Гиббса-Дюгема может быть реализовано с помощью концепции избыточной энергии Гиббса, т. е. превышения наблюдаемого уровня энергии Гиббса для смеси по отношению к величине, характерной для идеального раствора при тех же значениях температуры, давления и состава. Полная избыточная энергия Гиббса GE определяется соотношением 32 N G E  RT   ni  ln  i , (2.22) i 1 где R – универсальная газовая постоянная, Дж/(моль·К); Т – температура, К; ni – число молей i-ого компонента. Ключевой проблемой в расчете многокомпонентного фазового равновесия служит выражение для мольной энергии Гиббса, которое представляло бы собой аппроксимацию свойств смеси. При наличии выражения для избыточной энергии Гиббса можно получить выражения для расчета коэффициентов активности. Наибольшее распространение получили такие выражения для коэффициентов активности, как уравнение Вильсона, уравнение НРТЛ, уравнение ЮНИКВАК, метод ЮНИФАК и т. д. [12]. Рассмотрим в качестве примера задачу определения состава пара и температуры смеси по известному составу жидкости и давлению. Математическое описание включает в себя систему уравнений равновесия yi  P  xi   i  xi ,T   Pi 0 T  ; i  1, 2,3...N , (2.23) – парциальное давление компонента i и стехиометрические соотгде ношения xi i Pi 0 .  yi  1    P i 1 i 1 N N (2.24) В уравнениях (2.23) и (2.24) паровая фаза предполагается идеальной. Зависимость коэффициента активности в жидкой фазе γi от состава и температуры определяют одним из уравнений парожидкостного равновесия (Вильсона, ЮНИФАК и т. д.). Искомыми переменными в данном случае являются yi (N значений) и T (всего N + 1). Решение системы (2.23) и (2.24) проводится итерационными методами. Рассмотрим решение системы с применением алгоритма метода Ньютона [12]: 1. Задают начальное приближение температуры кипения T0. 2. Для заданного состава и температуры рассчитывают коэффициенты активности компонентов γi. 3. По уравнению (2.23) рассчитывают состав пара yi. 4. По уравнению (2.24) рассчитывают температуру кипения смеси Т. 5. Проверяют условие T  T '   . 6. Если условие выполняется, то расчет заканчивают. В противном случае расчет продолжается с п. 2. 33 Равновесие в системе «жидкость-жидкость» [8]. Описание равновесия между жидкими фазами во многом аналогично описанию системы «жидкость-пар». При расчете равновесия в системе «жидкость-жидкость» можно выделить две основные задачи:  определение состава равновесных фаз при заданной температуре по общему содержанию каждого компонента, присутствующего в смеси;  определение состава одной из равновесных фаз по заданному составу другой при известной температуре. Рассмотрим решение первой задачи. В качестве исходной информации для расчета состава равновесных фаз используются температура системы и количество молей каждого компонента. Неизвестными являются 2N составов, число молей в первой фазе М(1) и во второй фазах М(2), т. е. всего имеется 2N + 2 неизвестных. Математическое описание состоит из N уравнений материального баланса, N уравнений фазового равновесия и двух условий равенства мольных долей единице. В эквивалентной форме оно выглядит следующим образом:  (1) M F  zi  x ; i  1, 2,3...N ;  i (1) F     M 1 k M k   i i   (2) (1) (2.25)  xi  ki  xi ; i  1,2,3...N ;  N N  D   xi(1)   xi(2)  0.  i 1 i 1 Здесь ki – константа термодинамического равновесия между фазами; MF – общее количество молей в смеси; zi – общий состав смеси. Для совместного решения уравнений баланса и равновесия проводится следующая подстановка: xi( I )  xi(1) * xI( I )  x1( I ) *  i , i  2,3,..., N  1 , (2.26) где i  xi(2)  xi(1) * xI(2)  x1(1) * , i  2,3,..., N  1 . * Величины xi( j ) соответствуют концентрациям фаз, полученным в результате решения системы уравнения материального баланса при заданных константах равновесия. После такой подстановки решение проводится относительно следующего набора неизвестных: 34 x1  x1( I ) ; x2  x1(2) ; ................; xN  xN(2)1. Для решения системы нелинейных уравнений применим метод Ньютона. Алгоритм для расчета составов равновесных жидких фаз: 1. Задается число компонентов, температура, общий состав смеси zi и начальные оценки констант фазового равновесия ki. 2. Методом касательных рассчитывается число молей в первой фазе M(1). 3. Рассчитываются приращения по искомым переменным Δx1, Δx1, …, ΔxN. 4. Определяются новые значения искомых переменных x1  x1    x1; x2  x2    x1; ........................; xN  xN    xN . 5. 6. Оцениваются новые приближения констант фазового равновесия ki (2.25). Рассчитываются величины невязок rin системы уравнений фазового равновесия. Если выполняется условие ri  1; j , i  1,..., N , x j   2 , то расчет заканчивают и выводят результаты (xi(1), xi(2), M(1), M(2)). В противном случае продолжают расчет с п. 2 [8]. 2.3.2. Моделирование процесса массопередачи Детерминированное описание процесса переноса вещества в процессах массопередачи основано на фундаментальных законах диффузии Фика [7]. Для широко распространенных в промышленности процессов разделения, таких как процессы абсорбции, ректификации, экстракции, т. е. для процессов с так называемой свободной поверхностью раздела, существенно изменяющейся от взаимодействия двухфазных потоков, использование зависимостей, характеризующих детерминированные параметры, не приводит к желаемым результатам, и необходимо 35 прибегать к математическому моделированию, чтобы учесть стохастические (вероятностные) составляющие процессов [8]. В силу стохастического характера явлений массопереноса достижение равновесного состояния подчинено вероятностным законам распределения энергии и массы в пространстве и во времени. Главными причинами неравновесности в промышленных процессах являются: неравномерность распределения частиц потока по времени пребывания (по причине неравномерности профиля скоростей, турбулизации потоков, градиентов температуры и давления), обратный заброс фаз в результате механического уноса, недостаточное время контакта фаз. Поэтому при заданных конструктивных характеристиках аппарата время контакта фаз, определяемое гидродинамической структурой потоков, может оказаться недостаточным для достижения равновесия. Таким образом, важнейшим «элементарным» процессом при моделировании массообмена является процесс массопередачи. Рассмотрим основные уравнения массоотдачи и массопередачи. При отсутствии равновесия между фазами происходит перенос вещества из одной фазы в другую. Этот процесс называется массопередачей, которая является сложным процессом, состоящим из процессов переноса вещества в пределах каждой из фаз (массоотдача) и переноса вещества через границу раздела фаз [8]. Количество компонента i, переносимого в единицу времени t через поверхность F в единицу времени (закон Фика) составляет dC Wi   D  F  i . (2.27) dt В приведенном выражении D – коэффициент диффузии, м2/с; F – площадь поверхности массопередачи, м2; Ci – концентрация компонента i, моль. При рассмотрении уравнения массопередачи за движущую силу принимают разность между фактической концентрацией компонента в одной из фаз и равновесной концентрацией в ней данного компонента. Уравнение массоотдачи записывается в виде (2.28) Wi    F  , где Wi – количество вещества, переносимого в единицу времени; ∆ – движущая сила; β – коэффициент массоотдачи, представляющий собой количество вещества, переносимое внутри фазы в единицу времени через единицу поверхности при движущей силе, равной единице. В случае передачи вещества из паровой фазы с концентрацией x в жидкую y, уравнения массоотдачи запишутся в виде Wi   y  F  yi  yPi ; (2.29) Wi   x  F  xPi  xi ,   36   где , – коэффициенты массоотдачи в жидкой и паровой фазах соответственно; , – значения концентраций компонента i на поверхности раздела фаз в жидкой и паровой фазах соответственно. Исходя из условия равновесия фаз у поверхности их соприкосновения, после преобразований получим Wi  K y  F  ( yi  yi* ) , или Wi  K x  F  ( xi*  xi ), (2.30) где ∗ , ∗ – равновесные концентрации в жидкой и паровой фазах соответственно; Ky и Kx определяются уравнениями 1 1 m 1 1 1   ;   ; (2.31) K y  y  x K x m y  x и называются коэффициентами массопередачи, отнесенными к концентрации газа и жидкости соответственно. Уравнения (2.30) являются различными формами уравнения массопередачи [12]. 2.3.3. Моделирование процесса сепарации В химической технологии процесс сепарации (разделения) в основном газовой и жидкой фаз или жидких фаз, имеющих разные плотности, распространены очень широко [10, 17, 18]. Несмотря на большое разнообразие конструкций сепараторов, их можно условно разделить на два класса в соответствии с физическими принципами разделения газожидкостных смесей: гравитационные и инерционные [10]. На рис. 2.5 приведены примеры сепараторов. а б Рис. 2.5. Виды сепараторов: а – газосепаратор сетчатый; б – нефтегазовый сепаратор В гравитационных сепараторах, представляющих собой большие горизонтальные или вертикальные емкости, разделение фаз происходит за счет силы тяжести. Поскольку размеры капель, попадающих в сепа37 ратор из подводящего трубопровода, малы, то для их эффективного удаления из потока только за счет силы тяжести требуется длительное время, и, как следствие этого, сепараторы имеют большие размеры. В инерционных сепараторах разделение фаз происходит за счет сил инерции при обтекании газожидкостной смесью различных препятствий (сеток, струн и т. п.) и при закручивании потока в центробежных патрубках (циклонах). В современных конструкциях сепараторов используются оба принципа. Степень разделения газожидкостной смеси в сепараторах зависит от расхода газа, термобарических условий, а также от среднего радиуса капель, вносимых в сепаратор с потоком газа из подводящего трубопровода, который, в свою очередь, зависит от параметров трубопровода, а также от наличия установки предварительной конденсации перед сепаратором. Расчет процесса сепарации (однократного испарения) Считаем, что в процессе сепарации (рис. 2.6):  достигается состояние равновесия;  происходит однократное испарение компонентов смеси. Газ (G) Исходная cмесь (F) Жидкость (L) Рис. 2.6. Схема сепаратора Уравнение материального баланса процесса однократного испарения для многокомпонентной системы в целом можно представить как [17] F  G  L, (2.32) где F – количество исходного сырья, кг/ч; G – количество паровой фазы, кг/ч; L – количество жидкой фазы, кг/ч. Для i-го компонента системы материальный баланс запишется следующим образом: F  zi  G  yi  L  xi , 38 (2.33) где zi, xi, yi – мольные доли i-го компонента в исходном сырье и полученных жидкой и паровой фазах соответственно. В условиях равновесия yi  K i  xi ; (2.34) K i  Pi 0 / p, где Кi – константа фазового равновесия i-го компонента; Pi 0 – давление насыщенных паров i-го вещества, Па; P – общее давление, Па. Основное уравнение для расчета частичного однократного испарения многокомпонентной системы – zi  xi   , (2.35) 1  e  ( Ki  1) G – молярная доля пара (доля отгона) в конце процесса одноF кратного испарения. Контролем правильности решения уравнения (2.35) является выполнение условий: (2.36)  xi   yi  1. где e  Определить давление насыщенных паров компонентов можно по различным расчетным формулам, например: Антуана, Ашворта и др [15, 17]. В частности, формула Ашворта имеет следующий вид: Pi 0  105  exp 6,172  1  F T  F Ti    , (2.37) где Pi 0 – давление насыщенных паров, Па; T – температура однократного испарения, °C; Ti – температура кипения углеводорода, или средняя температура кипения углеводородной фракции, °C. Функцию F(T) находят из уравнения 1250 (2.38) F T    1. 2 T  273  108000  307,6 Подставляя в (2.38) Ti вместо Т, рассчитывают функцию F(Ti). Уравнение Антуана [15] – Bi ln Pi 0  Ai  , (2.39) T  Ci где Аi, Bi, Ci – коэффициенты уравнения; Т – температура процесса, Pi 0 – парциальное давление i-го компонента в системе. 39 2.3.4. Моделирование процесса ректификации Процесс ректификации – один из наиболее распространенных процессов разделения смесей в химической технологии. По определению, процессом ректификации называется термический способ разделения смесей путем многократного испарения и конденсации смеси, сопровождающийся тепло- и массообменом [7, 12]. Рассмотрим математическое описание ректификации. Модель тарельчатой ректификационной колонны основывается на следующих допущениях:  паровая фаза принимается идеальной;  жидкость на тарелках полностью перемешана;  количество тарелок – N;  смесь состоит из M компонентов;  исходное питание в количестве F состава zi подается на n-ю тарелку колонны. Сверху колонны отбирается дистиллят в количестве D состава xi,d, а снизу колонны кубовый продукт в количестве W состава xi,w. Схема ректификационной установки приведена на рис. 2.7. Рис. 2.7. Схема ректификационной колонны: Fn, D, W – потоки питания колонны, дистиллята и кубового остатка соответственно, кг/ч; L – поток флегмы, кг/ч; zi,n – состав питания 40 Схема потоков на тарелках колонны изображена на рис. 2.7,а. Рис. 2.7,а. Схема потоков жидкости и пара на тарелках ректификационной колонны Математическое описание включает следующие уравнения:  общего материального баланса на тарелках колонны: Gn 1  Ln 1  Fn  Gn  Ln  0, (2.40) где G, L – расходы пара и жидкости в колонне соответственно, кг/ч; n – номер тарелки колонны;  покомпонентного материального баланса: (2.41) Gn 1  yi ,n 1  Ln 1  xi ,n 1  Fn  zi ,n  Gn  yn  Ln  xn  0, где xi ,n , yi ,n – концентрации компонента i на тарелке n в жидкой и паровой фазах соответственно;  теплового баланса: Gn 1  H n 1  Ln 1  hn 1  Fn  hnF  Gn  H n  Ln  hn  0, (2.42) где Hn и hn – энтальпии парового и жидкостного потоков соответственно, Дж/моль; hnF – энтальпия потока питания, Дж/моль. Энтальпии парового H и жидкостного h потоков на каждой ступени разделения n определяются выражениями M H   yi  H i0 ; i 1 M h   xi  hi0 . (2.43) i 1 j  1,, M . Здесь H i0 , hi0 – стандартная энтальпия образования вещества j в паровой и жидкой фазах соответственно, Дж/моль; 41  фазового равновесия: yi*,n  Ki ,n  xi ,n , (2.44) где yi*,n – концентрация компонента в паре, находящегося в равновесии с жидкостью состава xi,n; Ki,n – константа равновесия между жидкостью и паром;  стехиометрические соотношения: M M  xn,i  1;  yn,i  1. i 1 (2.45) i 1 Если на тарелках колонны равновесие не достигается, то состав пара с предыдущей тарелки связывается с составом пара последующей через эффективность тарелки Ei n : Ei n  y n ,i  y n1,i y*n,i  y n1,i , (2.46) где y n,i , y n1,i – средняя концентрация компонента в потоках пара; y*n,i – концентрация компонента в паре, находящегося в равновесии с жидкостью состава xn. При расчете равновесия «жидкость-пар» отклонение от идеальности жидкой фазы учитывается с помощью коэффициента активности γ, определяемого как функция состава от температуры, например, по уравнению Вильсона [8]. С учетом уравнения (2.44) константу равновесия Ki,n можно записать в виде  i  xi  Pi 0 T  yi  . (2.47) P Зависимость от температуры давления насыщенных паров Pi0(T) i-го компонента определяется из уравнения [16] A (2.48) ln Pi 0 T   A1  2  A3  T  A4  ln(T ). T Система уравнений (2.39) – (2.48) нелинейная, и для ее решения должны применяться итерационные методы. Обычно для данной задачи различают два типа внешних условий [8]: 1. Внешние условия для решения задачи в проверочной постановке (расчет режимов работы колонны заданной конструкции):  состав и количество питания;  конструктивные параметры (диаметр колонны, число тарелок, межтарельчатое расстояние). 42 В результате решения задачи определяют: оптимальное флегмовое число, место ввода питания, состав продуктов разделения, профиль температуры по колонне. 2. Внешние условия для решения проектной задачи:  количество и состав разделяемой смеси;  содержание примесей в целевом продукте. Проектная задача является более общей и включает в себя проверочную. В результате решения задачи получают число тарелок в колонне, тарелку ввода питания, флегмовое число, диаметр колонны, межтарельчатое расстояние, тип тарелок, расход пара и жидкости, нагрузку на кипятильник и дефлегматор, состав продуктов разделения. Диаметр ректификационной колонны определяется из уравнения расхода [16] Dk  4V0 ,  w0 (2.49) где Dk – внутренний диаметр колонны; V0 – объемный расход пара в колонне; – допустимая скорость пара. Допустимая скорость пара в свободном сечении колонны рассчитывается таким образом, чтобы минимизировать унос флегмы паровым потоком на вышерасположенную тарелку: 1   L  G  2 4 (2.50) 0  8, 47  10 С1   ,  G  где L – плотность жидкости, кг/м3; G – плотность пара, кг/м3; С1 – эмпирический коэффициент, зависящий от расстояния между тарелками и поверхностного натяжения жидкости. По результатам расчета диаметра колонны выбирается межтарельчатое расстояние и тип тарелок. По данным составов дистиллята и кубового продукта, а также потоков пара с верха колонны и кубового отбора определяются тепловые нагрузки на кипятильник Qb и дефлегматор колонны Qc: M 0  Qb    rисп ,i  xWi   W ;  i 1  (2.51) M 0  Qc    rконд ,i  yNi   GN ,  i 1  где исп, – теплота испарения i-го компонента, Дж/моль; конд, – теплота конденсации i-го компонента, Дж/моль; – состав кубового остатка 43 колонны; – состав пара, уходящего с верхней тарелки колонны; GN – расход пара, уходящего с верхней тарелки колонны, кг/ч. Исходя из тепловых нагрузок кипятильника и дефлегматора, рассчитывают требуемые расходы греющего пара и охлаждающей воды. Для определения тарелки ввода питания можно использовать уравнение Керкбрайда [6] 0,206 2 N1 W  xiF xkW     , (2.52)   N 2  D  xiD xkF     где N1, N2 – число тарелок укрепляющей и исчерпывающей частей колонны соответственно; i – индекс тяжелого ключевого компонента, k – индекс легкого ключевого компонента. 2.3.5. Моделирование процесса абсорбции Абсорбцией называют процесс поглощения газа или пара жидким поглотителем абсорбентом. Рассмотрим процесс абсорбции в насадочной колонне (рис. 2.8.) [19]. y2 Ll x1, Gg, y1 1 x2 Рис. 2.8. Насадочная абсорбционная колонна Будем считать температуру по высоте колонны постоянной и равной заданному значению. Вследствие этого из уравнений математического описания можно исключить уравнения тепловых балансов. Содержание абсорбируемых компонентов в газе относительно мало, и допускается, что условия фазового равновесия могут быть описаны для всех компонентов с использованием закона Генри [9]: y   K i  xi , где Ki – постоянная Генри для каждого компонента; yi. – равновесная концентрация i-го компонента в паровой фазе; – концентрация i-го компонента в жидкой фазе. 44 Материальный баланс массообменного процесса абсорбции можно представить как разность потоков на входе и выходе аппарата для паровой и жидкой фаз: Gg   y1  y2   Ll   x2  x1  , (2.53) где Gg – расход инертного газа, м3/с; Ll – расход жидкости, м3/с. При абсорбции процесс массопередачи происходит на поверхности соприкосновения фаз, поэтому в абсорберах необходимо обеспечить наибольшую поверхность их соприкосновения. В зависимости от конструктивных особенностей абсорберов процессы, протекающие в этих аппаратах можно описать различными гидродинамическими моделями. Если в аппаратах колонного типа потоки газа и жидкости движутся в режиме идеального вытеснения, то в стационарном режиме модель процесса абсорбции будет представлять собой следующую систему уравнений материального баланса:  для газовой фазы dy    g S y  y* ; Gg (2.54) dl  для жидкой фазы dx  l  x*  x , Gl (2.55) dl где  g ,  l – коэффициенты массоотдачи для газовой и жидкой фаз, с–1; S – поверхность массообмена, м2; L – высота аппарата, м. y(0) = y0; при l = L x(l) = x0. Граничные условия: при l = 0 Если структура потоков в аппарате не соответствует режиму идеального вытеснения, то могут быть применены другие модели, например ячеечная модель.     2.3.6. Моделирование процесса адсорбции Адсорбция – процесс поглощения газов, паров или жидкостей поверхностью пористых твердых тел. Процессы адсорбции широко применяются для очистки и осушки газов, разделения смесей (газов и паров), регенерации растворителей, очистки от примесей [19]. Для проведения данного процесса применяют адсорберы следующих типов:  c неподвижным зернистым адсорбентом;  с движущимся слоем адсорбента;  с кипящим (псевдоожиженным) слоем. 45 Из-за накопления сорбата на поверхности сорбента свойства последнего постоянно изменяются, и, следовательно, весь процесс адсорбции в целом является нестационарным. Так как концентрация адсорбирующегося вещества меняется по высоте слоя сорбента, уравнения материального и теплового баланса можно записать лишь для элементарного объема. В результате математическое описание представляет собой систему дифференциальных уравнений в частных производных. На основании исследований процесса адсорбции было установлено, что для подвижной фазы можно воспользоваться диффузионной моделью с учетом продольного перемешивания, кроме того, в подвижной фазе происходит массообмен между фазами. Математическое описание для неподвижной фазы включает составляющую, описывающую массообмен между фазами. Тогда двухфазную математическую модель процесса адсорбции можно представить следующим образом:  уравнение материального баланса для подвижной фазы – C  2C C  D 2 u    C  C* ; t l l (2.56)  уравнение материального баланса для неподвижной фазы – dCн    C  C*  , (2.57) dt где С, Сн – концентрации компонентов в подвижной и неподвижной фазах. Начальные условия: при t = 0 C(0,l) = C0, Cн = 0. C  0. Граничные условия: при l = 0 C = C0, l Вопросы для самоконтроля 1. 2. 3. 4. 5. 6. Назовите основные массообменные процессы, применяющиеся в химической технологии. Какие фундаментальные законы лежат в основе описания массообменных процессов? Что такое фазовое равновесие? Какие методы расчета констант фазового равновесия вы знаете? Какие основные задачи решаются при моделировании равновесия «жидкость-пар»? Как выражается условие термодинамического равновесия между жидкостью и паром? В системе «жидкость-жидкость»? Какие вы знаете соотношения, связывающие активность компонента с составом смеси и температурой? 46 7. 8. 9. 10. 11. 12. 13. 14. 15. Что такое массопередача и массоотдача? Как связаны между собой коэффициенты массоотдачи и массопередачи? Что такое ректификация? Какие уравнения входят в математическое описание процесса ректификации? Что является исходными данными и результатом расчета при моделировании процесса ректификации? В чем коренное отличие моделирования насадочной колонны от тарельчатой? Какие численные методы, применяющиеся для решения систем нелинейных уравнений, вы знаете? В чем заключается различие процессов сепарации и ректификации? Какими математическими моделями описывается процесс абсорбции? Какими математическими моделями описывается процесс адсорбции? 2.4. Математическое моделирование кинетики химических реакций 2.4.1. Основные понятия химической кинетики Учение о скоростях химических реакций называется химической кинетикой. Химическая кинетика как наука начала формироваться в 50–70 гг. XIX в. Скорость химической реакции есть изменение числа молей реагентов в результате химического взаимодействия в единицу времени в единице объема (для гомогенных реакций) или на единице поверхности (для гетерогенных процессов) [20]. В соответствии с этим скорость реакции записывается следующим образом: 1 dN W=  ; для гомогенной реакции (2.58) V dt 1 d C V  W   ; V dt 1 dN W=  , для гетерогенной реакции (2.59) S dt где V – объем реакционной фазы, м3; N – количество молей вещества; С – концентрация реагента, моль/м3; t – время, с.; S – поверхность катализатора, м2. Выражение (2.58) может быть преобразовано к виду  VdC C dV   W    Vdt V dt 47  .  (2.60) Для реакций, идущих при постоянном объеме, второе слагаемое в уравнении (2.60) равно нулю, тогда скорость реакции dC W  . (2.61) dt Для реакторов периодического действия, в которых концентрации реагирующих веществ в каждой точке реакционного объема в ходе реакции непрерывно изменяются во времени, скорость химической реакции есть количество молей данного вещества, реагирующее в единицу времени в единице объема: dN 1 Wi   i  (2.62) dt V или на единицу поверхности для гетерогенных каталитических реакций: dN 1 dN 1 Wi   i    i  , (2.63) dt S dt  0V где Ni – текущее количество i-го компонента реакционной смеси, моль; V – объем реакционной смеси или слоя катализатора (объем реактора), м3; 0 – удельная поверхность катализатора, м2/м3. Для реакторов непрерывного действия полного вытеснения, в которых при установившемся режиме концентрация вещества непрерывно изменяется по длине аппарата, скорость химической реакции есть количество молей проходящего через реактор в единицу времени вещества, реагирующего в единице объема [1, 4, 21]: dn dni , Wi   i   (2.64) dV d  v где ni – мольный расход i-го компонента реакционной смеси, моль/с; v – объемная скорость подачи реакционной смеси, м3/с;  – время контакта, с. Для реактора непрерывного действия полного смешения, при установившемся режиме, n n Wi  i 0 i , (2.65) V где ni0, ni – начальное и конечное количество i-го компонента реакционной смеси, моль/с. На практике обычно измеряют скорость изменения мольной концентрации Сi (моль/м3):  для реактора периодического действия: n Ci  i ; V dni  Ci dV  VdCi . 48  для реактора непрерывного действия: Ci  ni ; v dni  Ci dv  vdCi . Если реакция не сопровождается изменением объема, то для реактора идеального вытеснения vdCi dCi dC dx    i  Ci 0 i . Wi   (2.66) dV d V / v  d d Для реактора непрерывного действия идеального смешения C v  Ci v Ci 0  Ci Ci 0  Ci x Wi  i 0    Ci 0 i , V V v   где xi – степень превращения: xi  (2.67) ni 0  ni ;  – среднее время пребывания, с; ni 0  = V/v. В 1862–1867 гг. норвежские ученые Гульдберг и Вааге дали начальную формулировку закона действующих масс. Согласно кинетическому закону действующих масс скорость элементарной реакции при заданной температуре пропорциональна произведению концентраций реагирующих веществ в степенях, равных их стехиометрическим коэффициентам [3, 4, 22]. Тогда при протекании химической реакции k    A B c D по закону действующих масс выражение скорости запишется следующим образом: A W  k С   CBB , (2.68) A где Ci – концентрация i-го вещества, k – константа скорости; υA, υB – стехиометрические коэффициенты веществ А и В. Уравнение (2.68) справедливо для элементарных реакций. Между скоростями реакции по отдельным компонентам (обозначим их WA, WB, …) и общей скоростью реакции W существует следующее стехиометрическое соотношение: W W W W (2.69) W  A  B  C  D.   A B c D 49 Чтобы применить закон действующих масс к сложной химической реакции, необходимо представить ее в виде элементарных стадий и применить этот закон к каждой стадии отдельно. Химическая кинетика в полной мере была сформулирована в работах Вант-Гоффа и Аррениуса в 80-х гг. XIX в.; был разъяснен смысл порядков реакций, введено понятие энергия активации, моно-, би- и полимолекулярных реакций. Вант-Гофф, а впоследствии Аррениус, развивший его идеи, утверждали, что «температура не есть причина реакции, температура – причина изменения скорости реакции» [20]: k  A  e  E RT , (2.70) где А – предэкспоненциальный множитель; Е – энергия активации, Дж/моль; R – газовая постоянная, Дж/моль·К; Т – температура, К. При заданных внешних условиях (температуре, давлении, составе, среде, в которой протекает реакция) скорость реакции является функцией концентраций реагирующих веществ: (2.71) W  kf (C А , СС ,..., СZ ). Сопоставляя уравнения (2.61) и (2.71), получаем кинетическое уравнение dCA dC A  kf .  kf (C А , СС ,..., СZ ) (2.72) dt dt Уравнение, отражающее изменение концентрации вещества во времени в ходе химической реакции, называется кинетическим, а кривая c  f  t  – кинетической кривой. 2.4.2. Моделирование кинетики гомогенных химических реакций Кинетические уравнения связывают скорость реакции с параметрами, от которых она зависит, например: с концентрацией, температурой, давлением, активностью катализатора. Исследование кинетических закономерностей протекания химических реакций методом математического моделирования заключается в определении изменения концентраций реагирующих веществ во времени. Пусть протекают химические реакции k1 k2 A   2C  D. На основании закона действующих масс запишем уравнения скоростей химических реакций и составим кинетическую модель: 50 W1  k1  C A ; W2  k 2  CC2 ; dC A =-k1C A ; (2.73) dt dCD  k2 CBCC2 ; dt dCC =2k1CA  2k2CC2 , dt где СА, СС, СD – концентрации веществ, моль/л. Систему обыкновенных дифференциальных уравнений (ОДУ) первого порядка (2.73) можно решить с использованием алгоритмов численных методов решения ОДУ (Эйлера, Рунге-Кутты [11]). На рис. 2.9, 2.10 приведены результаты исследования влияния температуры на степень превращения исходного реагента и на концентрацию веществ с применением математической модели (2.73). Полученные результаты позволяют сделать вывод об оптимальном времени проведения процесса с целью получения целевого продукта. Математическая модель (2.73) также позволяет исследовать влияние состава сырья на выход продуктов реакции. Известно, что скорость химической реакции зависит от температуры, поэтому, чтобы исследовать влияние температуры на ход процесса, необходимо в кинетической модели (2.73) константу скорости химической реакции выразить по уравнению Аррениуса (2.70). C, моль/л 1 0,8 0,6 0,4 0,2 t, с 2 Са 4 Сс 6 Сд 8 10 Рис. 2.9. Изменение концентрации реагирующих веществ 51 с течением времени Степень превращения,% 100 95 90 85 80 75 550 560 570 580 590 600 610 Температура, оС Рис. 2.10. Зависимость степени превращения исходного вещества от температуры 2.4.3. Моделирование кинетики гетерогенных химических реакций Основы гетерогенной химической кинетики заложены в работах Лэнгмюра, Темкина и др. [20, 21, 23, 24]. В этих работах сформулировано понятие идеального адсорбированного слоя, базирующееся на аналогии с представлениями гомогенной кинетики. Эта модель использует следующие предположения:  равноценность всех участков поверхности катализатора и независимость энергии хемосорбции от степени заполнения поверхности различными адсорбентами;  неизменность катализатора и независимость его свойств от состава реакционной смеси и ее воздействия на катализатор;  равновесное распределение энергии. Формальным аналогом кинетического закона действующих масс для элементарных процессов на твердых поверхностях является закон действующих поверхностей (ЗДП) [20, 21]. Согласно его первоначальной формулировке скорость химической реакции пропорциональна произведению поверхностных концентраций реагирующих веществ в степенях, равных стехиометрическим соотношениям, в которых они вступают во взаимодействие (2.74): 1 2  Z0 W  kQZ1  QZ2  ...  Q i 2  W  k   z1   z2  ...   52  z0 n  k  zii   z0 . (2.74) i 1 где k – константа скорости;  Zi – доля поверхности, занятой i-й адсорбированной частицей;  Z0 – доля свободной поверхности; i – стехиометрические коэффициенты стадий;  – изменение числа молей при протекании химической реакции. Пусть протекает элементарная химическая реакция 1 A1   2 A2  ...  1' A1'   2' A2'  ... . (2.74а) При этом все вещества вступают во взаимодействие из адсорбированного состояния. Обозначим: zi – доля поверхности, занятая i-м адсорбированным веществом. Тогда, в соответствии с законом действующих поверхностей, скорость необратимой реакции (2.74а) можно записать как n W  k  zii  k   z1i   z22   z0 . (2.75) i 1 Если не все вещества вступают во взаимодействие из адсорбированного состояния, а реагируют непосредственно из газовой фазы, то в более общем виде выражение закона действующих поверхностей можно представить следующим образом: n m W  k   zi  C j j , i i 1  (2.76) j 1 где C j – парциальные давления (концентрации) j-веществ, реагирующих из газовой фазы; n, m – количество веществ, адсорбированных на поверхности катализатора и реагирующих из газовой фазы. Например, пусть протекает адсорбция водорода на активном центре катализатора Z с образованием адсорбированного поверхностного соединения ZH2: k H 2  Z   ZH 2 , тогда на основании ЗДП скорость данной элементарной химической реакции можно записать как W  k Z  C . H2 В качестве основного фактора, определяющего кинетические зависимости, вначале рассматривался фактор вытеснения, «борьбы» компонентов реакционной смеси за места на поверхности катализатора. При этом принималось дополнительное предположение о высокой скорости адсорбционных и десорбционных стадий по сравнению с собственно химическими превращениями. 53 Последующие исследования показали существенную ограниченность этих предположений. Тем не менее Хиншельвудом, Швабом, Хоугеном, Ватсоном и другими получены уравнения, удовлетворительно описывающие кинетический эксперимент в определенном интервале изменения параметров [24]. Общая формула кинетического уравнения, соответствующего этим предположениям, имеет следующий вид: vi n n     W   k   Ci  1   k pi  Ci  , i 1   i 1   где k – константа скорости; Сi – концентрация i-го реагента газовой среды; k pi – константа равновесия стадии адсорбции i-го компонента; i – стехиометрические коэффициенты. Пример Рассмотрим сложную гетерогенную химическую реакцию гидрокрекинга толуола. Детальный механизм гетерогенной химической реакции: Здесь Z – активные центры на поверхности катализатора; ZH2 и т. д. – адсорбированные промежуточные соединения. Запишем скорости элементарных стадий механизма по закону действующих поверхностей: r1  k1  CH   Z ; 2 r1  k1   ZН ; 2 r2  k 2  CC7 H8   ZH 2 ; r3  k3   ZC7 H8 H 2 ; r3  k 3  CC6 H 6  CCH 4   Z . Математическая модель данного химического процесса будет представлять собой систему дифференциальных уравнений, выражающих изменение концентраций наблюдаемых веществ и промежуточных соединений во времени: 54  dCH2  dt  r1  r1;   dCC7H8  dt  r2 ;   dCCH4  dt  r3  r3 ;   dCC6H6  r3  r3 ; (2.77)  dt   d Z  dt  r1  r1  r3  r3 ;   d ZН2  dt  r1  r1  r2 ;   d ZC7H8 H 2  r2  r3  r3 .  dt  При решении системы дифференциальных уравнений (2.77) можно использовать численные методы Эйлера и Рунге-Кутты. Методы построения кинетических моделей гетерогенных химических реакций Метод, основанный на элементах теории Ленгмюра Сформулируем алгоритм построения кинетической модели гетерогенной химической реакции методом Лэнгмюра. 1. Записывается механизм гетерогенной химической реакции в виде совокупности элементарных стадий. 2. Записываются выражения скоростей реакций элементарных стадий на основании закона действующих поверхностей. 3. Определяется лимитирующая стадия. Все остальные стадии являются быстрыми и равновесными. 4. Для равновесных стадий приравниваются скорости в прямом и обратном направлениях и концентрации промежуточных соединений выражаются через концентрацию промежуточного соединения, участвующего в лимитирующей стадии и концентрации наблюдаемых веществ. 5. В результате решения системы алгебраических уравнений определяются выражения для концентрации промежуточного соединения, входящего в уравнение скорости лимитирующей стадии, через концентрации наблюдаемых веществ и константы скоростей эле55 ментарных стадий. При этом используется уравнение нормировки, в котором количество активных центров на поверхности катализатора принято равным единице и равно сумме долей поверхности катализатора, занятых промежуточными соединениями и свободными центрами. 6. Записывается выражение скорости гетерогенной химической реакции по лимитирующей стадии. Пример. Рассмотрим вывод выражения скорости для следующей гетерогенной химической реакции: По закону действующих поверхностей запишем скорости элементарных стадий: r1  k1   z  C A ; r1  k 1   zA ; r2  k 2   zA ; r3  k3   zB ; r3  k 3   z  C B . Примем вторую стадию за лимитирующую и запишем выражение скорости стадии по ЗДП: (2.78) W  r2  k 2   zA . В этом выражении необходимо выразить концентрацию промежуточного соединения θzA через концентрации наблюдаемых компонентов. Все остальные стадии будем считать быстрыми и равновесными. В этом случае скорость прямой реакции равна скорости обратной: k1   z  C A  k 1   zA ; (2.79) k3   zB  k 3  C B   z . Выразим  z : z  k1   zA . k1  C A 56 Запишем уравнение нормировки  zA   zB + z  1 и решим систему алгебраических уравнений (2.79) относительно θzA:   zA 1    zB  k1  k1k3CB     1; k1C A k1k3C A  k3CB z k1k3CB zA  ; k3 k1k3C A  zA  1 k k k C 1  1  1 3 B k1C A k1k3C A ;  zA  k1k3C A ; k1k3C A  k1k3  k1k3CB W k1k2 k3C A . k1k3C A  k1k3  k1k3CB Метод стационарных концентраций Стационарным считается такой режим, при котором параметры в каждой точке системы в заданный промежуток времени сохраняются постоянными. Так как скорость химического процесса существенно зависит от концентраций промежуточных веществ, то при стационарном режиме ее протекания концентрации промежуточных соединений должны оставаться постоянными. Таким образом, стационарность протекания химической реакции означает стационарность концентраций промежуточных соединений, которые должны успевать многократно возникнуть и разложиться в ходе химической реакции, уступая место другим промежуточным соединениям. Тогда условие стационарности можно выразить следующим образом [23]: di  0. dt Выполнение условия стационарности означает, что алгебраическая сумма скоростей образования и расходования каждого из промежуточных соединений равна нулю. Для реализации стационарного режима необходимо сохранение постоянства внешних параметров режима, а для гетерогенных каталитических реакций – и постоянство каталитической активности катализатора. Сформулируем алгоритм построения кинетической модели гетерогенной химической реакции методом стационарных концентраций: 57 1. Записывается механизм гетерогенной химической реакции в виде совокупности элементарных стадий. 2. Записываются выражения скоростей реакций элементарных стадий на основании закона действующих поверхностей. 3. Записывается условие стационарности по промежуточным соединениям. 4. Решается полученная система алгебраических уравнений и определяются выражения для концентраций промежуточных соединений. 5. Записывается выражение скорости гетерогенной химической реакции по любой стадии. Пример. Составить выражение скорости для гетерогенной химической реакции. Механизм химической реакции: k1 1. A+z 2. zA k2 k -1 k -2 zA z+B А=В Скорости элементарных реакций: W1  r1  r-1 ; W2  r2  r2 . Запишем выражения для скоростей элементарных стадий по закону действующих поверхностей: r1  k1  C A   z , r2  k 2   zA ; r1  k 1   zA , r2  k 2  C B   z . Условие стационарности: dCz  0; dt dCzA  0; dt  dCz  dt  r1  r1  r2  r2  0;   dCzA  r  r  r  r  0. 1 2 1 2  dt Используя уравнение нормировки  z   zA  1, 58 (2.80) преобразуем систему уравнений (2.80) к следующему виду:  k1C A z  k 1 zA  k 2 zA  k 2C B z  0;  zA  k 1  k 2    z  k1C A  k 2C B  ; z   k1  k2  zA ; k1CA  k2CB  zA 1    zA   zA  r2  k1  k2    1; k1C A  k2CB  1 k1  k2 1 k1C A  k2CB ; k1C A  k2CB ; k1C A  k2CB  k1  k2 k2  k1C A  k 2  k2CB k1C A  k2CB  k1  k2 . Так как реакция обратимая, находим выражение скорости r–2: z  k1  k2 (k1C A  k2CB )  . (k1C A  k2CB ) k1C A  k2CB  k1  k2 В результате преобразований получим z  k1  k2 . k1C A  k2CB  k1  k2 Подставим полученное выражение в r–2: r2  k2  k1CB  k2  k2CB . k1C A  k2CB  k1  k2 В результате получим выражение скорости реакции W =r2 -r-2 = k2  k1C A  k2  k2CB k k C  k k C  2 1 B 2 2 B ; k1C A  k2CB  k1  k2 k1C A  k2CB  k1  k2 W= k1k2C A  k1  k2CB . k1C A  k2CB  k1  k2 Если в механизме есть хоть одна необратимая стадия, то весь процесс становится необратимым. 59 Метод построения кинетических моделей с использованием элементов теории графов Рассмотрим вывод выражения скорости на примере следующей реакции: . Механизм реакции представим в виде последовательности элементарных стадий: Граф для данного механизма реакции представим в виде Всякое графическое изображение конечного множества элементов и взаимосвязей между ними называется графом. Граф состоит из вершин (узлов) и ребер (или дуг)[20]. Дуга – это линия, соединяющая вершины графа. Цикл – замкнутая совокупность дуг. Дерево граф (каркас) – незамкнутая совокупность дуг, проходящих через все вершины графа и входящих в данную вершину. Например, дерево граф для вершины zC: 60 Вес дуги графа – это отношение скорости химической реакции данной элементарной стадии к концентрации промежуточного вещества, из которого эта дуга выходит: b1  b2  k1C A z z k2CB zA b3   zA k3 zC  zC  k1C A ;  k2CB ;  k3 . Вес дерева графа равен произведению весов дуг, его образующих: Bz  b2  b3  k 2 k3C B ; BzA  b1  b3  k1k3C A ; BzC  b1  b2  k1k 2C AC B . По определению, концентрация промежуточного вещества есть отношение веса каркаса для данного промежуточного вещества к сумме весов всех остальных каркасов: z  b2b3 k2 k3CB  b2b3  b1b3  b1b2 k2 k3CB  k1k3C A  k1k2C ACB . В основе данного метода лежит теория стационарного протекания химической реакции. Следовательно, скорость брутто-реакции можно записывать по любой из элементарных стадий механизма, например: W  r1  k1C A z ; W k1k2 k3CBC A k2 k3CB  k1k3C A  k1k2C ACB . В графе механизма линейной одномаршрутной реакции, содержащем n вершин, каждая вершина может иметь один прямой и один обратный каркасы и (n–2) смешанных каркаса. Общее количество каркасов Nk  2n  n  n  2  n2 . Пусть протекает обратимая химическая реакция 61 Граф механизма данной реакции будет иметь следующий вид: Веса каркасов: прямые обратные смешанные + + – – b2 b1 ; b3+ b1–; I b2 b3 ; I b1+ b3+; b3– b2–; b1+ b2–; I b1+ b2+. b1– b2–. b2+ b3–. В общем случае вес дуги графа определяется отношением скорости стадии к концентрации промежуточного вещества, которое в ней участвует: vs ; b  zs vs  bs  , zs 1  s где bs , bs – веса дуг в прямом и обратном направлениях соответственно;   vs , vs – скорости реакций в прямом и обратном направлениях; z s , z s 1 – концентрация промежуточного вещества в стадиях s или s + 1. В целом выражение скорости гетерогенной химической реакции можно записать следующим образом: W n n i 1 i 1  bi   bi B прi   Bобрi   Bсмешi . (2.81) Скорость химической реакции равна произведению весов дуг в прямом направлении «минус» произведение весов дуг в обратном направлении, отнесенное к сумме весов каркасов в прямом, обратном и смешанном направлениях. Вопросы для самоконтроля 1. 2. Каковы основные концепции формальной кинетики? Что такое скорость химической реакции? Как она определяется? 62 3. 4. Какова температурная зависимость скорости химической реакции? Какой закон лежит в основе формальной кинетики? Дайте его формулировку. 5. Какие численные методы используются для решения кинетических уравнений? 6. Какие экспериментальные данные необходимы для оценки кинетических констант и энергий активации? 7. Какова физическая природа многостадийного протекания гетерогенной химической реакции? 8. Перечислите особенности применения закона действующих поверхностей и назовите его отличие от закона действующих масс. 9. Перечислите методы построения кинетических моделей гетерогенных химических реакций. 10. Дать характеристику методу, основанный на использовании изотермы Ленгмюра. 11. Дать характеристику методу стационарных концентраций. 12. Дать характеристику методу с использованием теории графов. 2.5. Моделирование гомогенных химических реакторов Одним из основных элементов любой химико-технологической системы (ХТС) является химический реактор. Химическим реактором называется аппарат, в котором осуществляются химические процессы, сочетающие химические реакции с массо- и теплопереносом, с целью получения определенного вещества. Типичные реакторы – это контактные аппараты, реакторы с механическим, пневматическим и струйным перемешиванием, промышленные печи и т. д. От правильности выбора реактора и его совершенства зависит эффективность всего технологического процесса [21]. В химическом реакторе имеют место большое количество различных явлений и их взаимодействия. Одновременный их учет может привести к усложнению математической модели и трудностям, связанным с ее использованием [5]. Поэтому при разработке математических моделей химических реакторов используется системный подход, который реализуется в иерархической схеме построения модели реактора. Проводится декомпозиция процесса на составляющие. Процесс на более низком масштабном уровне является одной из составляющих более высокого уровня. Поэтому математическая модель процесса в целом представляет собой синтез моделей явлений разного масштаба. Реальные химические реакторы существенно отличаются друг от друга, следовательно, задача построения математических моделей должна решаться в каждом конкретном случае с учетом особенностей процесса и его конструктивного оформления. 63 2.5.1. Классификация реакторов В химической технологии применяют всевозможные типы реакторов, имеющие существенные различия [21]. Тем не менее реакторы можно классифицировать по некоторым признакам: 1. В зависимости от фазового состояния реагирующих веществ реакторы могут быть гомогенными или гетерогенными. 2. По характеру операций загрузки и выгрузки различают реакторы периодического, непрерывного и полупериодического действия. 3. По режиму движения реакционной среды или по структуре потоков вещества:  реакторы идеального смешения;  реакторы идеального вытеснения;  реакторы с продольным перемешиванием;  реакторы с продольным и радиальным перемешиванием;  реакторы с комбинированной структурой потока. 4. По тепловому режиму реакторы разделяются на изотермические, адиабатические и политропические. Изотермические реакторы имеют одну постоянную температуру во всех точках реакционного пространства. Адиабатический реактор не имеет теплообмена с окружающей средой. Это достигается хорошей тепловой изоляцией. В политропическом реакторе происходит теплообмен с окружающей средой. 5. По конструктивным признакам различают емкостные, трубчатые, комбинированные. Приведенная классификация свидетельствует о том, что реальные химические реакторы характеризуются большим числом свойств, поэтому при построении математической модели химического реактора необходимо выделить и учесть наиболее важные свойства, т. к. учесть одновременно все свойства невозможно. 2.5.2. Математическая модель реактора идеального смешения Математическое описание реактора идеального смешения (рис. 2.11,а и б) характеризует изменение концентраций в реакционной среде во времени, которое обусловлено движением потока (гидродинамический фактор) и химическим превращением (кинетический фактор). Поэтому модель реактора идеального смешения можно построить на основе типовой модели идеального смешения формула (2.2) с учетом скорости химической реакции [1–5]. С учетом кинетического фактора динамическая модель изотермического реактора идеального смешения непрерывного действия будет иметь вид 64 dСi 1    Cвх  Свых   Wi . (2.82) dt  где Сi – концентрация i-го вещества, кмоль/м3;Wi – скорость реакций по i-му веществу, кмоль/м3. Рис. 2.11,а. Реактор с мешалкой v, Свх V v, Свых Рис. 2.11,б. Схема реактора идеального смешения Уравнение (2.82) записывается по каждому из компонентов, участвующих в реакции. Система приведенных уравнений материального баланса (2.82) является математической моделью реактора идеального смешения. Запишем математическую модель реактора идеального смешения для реакции k A  B 65 dC A 1  C A0  C A  kC A ; dt  (2.83) dCB 1  CB0  CB  kC A . dt  Начальные условия: при t = 0 С A (0)  C A0 ; C B (0)  0. Это система уравнений материального баланса (2.83) для динамического режима работы реактора. В стационарном режиме работы аппарата dС A dCB  0;  0. dt dt При решении данных уравнений можно найти следующие основные параметры:  время контакта, характеризующее объем аппарата;  степень превращения и селективность процесса;  изменение концентраций реагирующих веществ как функцию от времени контакта; 1 (C A0  C A )  kC A  0;      C A0  C A   k C A ;  C A0  C A CA  xA  kC A CA0 ; 1 k C A0  C A C A0 ; ; C A  C A0 (1  xA );  C A0  C A0 (1  xA ) kC A0 (1  xA ) ; xA . k (1  xA ) Аналогично уравнению материального баланса реактора идеального смешения (2.82) записывается уравнение теплового баланса. Так, для адиабатического реактора получим N  смСpсм см см dT  Ср   Т вх  Т     H i   Wi , (2.84) dt  j 1  66 а для политропического реактора N Сpсм см dT Ср   Т вх  Т     H i   Wi  KF T , dt  j 1 (2.85) где Wi – скорость i-й химической реакции; ∆Hi – тепловой эффект i-й химической реакции; Сpсм – теплоемкость реакционной смеси; Т вх – температура на входе в реактор; Т – текущее значение температуры. Теплоемкость i-го вещества как функция температуры описывается следующим уравнением: С pi  (ai  bi  T  ci  T 2  di  T 2 )  4,1887. (2.86) Теплоемкость смеси вычисляется по правилу аддитивности: см p C N   C pi  Ci , (2.87) i 1 где Сi – концентрация i-го вещества смеси, м. д. Зависимость константы скорости химической реакции от температуры выражается уравнением Аррениуса (2.70). Для того чтобы исследовать работу реактора идеального смешения в динамическом режиме работы, т. е. проследить изменение концентрации реагирующих веществ и температуры во времени и на выходе из реактора, необходимо решить систему дифференциальных уравнений материального баланса по каждому из компонентов совместно с уравнением теплового баланса. 2.5.3. Математическая модель реактора идеального вытеснения Математическое описание реактора идеального вытеснения характеризует изменение концентраций в реакционной среде во времени, которое обусловлено движением потока (гидродинамический фактор) и химическим превращением (кинетический фактор). Поэтому модель реактора идеального вытеснения можно построить на основании типовой гидродинамической модели идеального вытеснения формула (2.4) с учетом скорости химической реакции [1–5]. С учетом кинетического фактора динамическая модель изотермического реактора идеального вытеснения непрерывного действия будет иметь вид Ci C  u i  Wi , (2.88) t l где Сi – концентрация соответствующего i-го вещества; Wi – скорость реакции по i-му веществу. 67 Уравнение теплового баланса адиабатического реактора идеального вытеснения: N см см T см см T (2.89)   Ср   U    Ср    H  W . i i t l i  1 Уравнение (2.89) записывается по каждому из компонентов, участk вующих в реакции. Например, для реакции А   B , протекающей в изотермическом реакторе идеального вытеснения, математическая модель (динамический режим) будет иметь вид C C A  u  A  k  C ; A t l (2.90) C C B  u  B  k  C . A t l В установившемся (стационарном) режиме работы реактора C A  0; t (2.91) С B  0, t тогда уравнения (2.90) примут следующий вид: dC u A  k  C A ; dl (2.92) dCе u  k  CA. dl l Так как   , то уравнения (2.92) примут вид u dC A  k  C A ; d (2.93) dCB  k  CA , d где  – время контакта, с. При решении системы уравнений (2.93) можно найти следующие параметры:  время контакта;  степень превращения и селективность процесса;  изменение концентраций реагирующих веществ как функцию от времени контакта;  68  C 1 A dC A    ; K CA C A   1 (nC A0 ); K   C 1 1 CA n A ;  n 0 ; K C A0 K CA CA  CA0 e K ; C A  C A0 (1  x A );  C A0 1 ; n K C A0 (1  x A )  1 1 . n K 1  xA Для того чтобы с применением данной модели выполнить исследования химико-технологического процесса, протекающего в реакторе идеального вытеснения, необходимо решить систему дифференциальных уравнений (2.93). 2.5.4. Исследование химического процесса, протекающего в гомогенном реакторе идеального смешения Пусть в реакторе идеального смешения протекает химическая реакция превращения н-октана в и-октан и в продукты крекинга: - H1  2 н  С8H18  изо  C8H18 C4H10  C4 H8 , k1 где H 1 = –7,03 Дж/моль при 700 К – экзотермическая реакция; H 2 = +85,89 Дж/моль – эндотермическая реакция. Представим химическую реакцию в виде k1 k2 A   B  C  D . Математическая модель процесса, с учетом уравнений (2.82, 2.84), может быть записана в виде следующей системы уравнений материального и теплового балансов: 69   dC A 1   C A0  C A  k1  C A ; dt  dCB 1  CB0  CB  k1  C A  k2  CB ; dt  dCC 1  CC0  CC  k2  CB ; dt  dCD 1  (CD0  CD )  k2  CB ; dt  (Q1  k1  C A  Q2  k2  CB )  R  T / p dT 1  T0  T   , dt  Cp     (2.94) где Р – давление в реакторе, Мпа. Начальные условия: при t = 0 CА (0) = CА,0;CB(0) = CC(0) = CD(0) = 0, R – универсальная газовая м3  МПа . постоянная, R  0,00845 кмоль  К Так как тепловой эффект реакции (Qi) равен величине энтальпии i-й реакции ( H i ) с обратным знаком: Qi  H i , то Qi  7,03 Дж/моль , Q2  85,89 Дж/моль . Для решения системы дифференциальных уравнений может быть использован численный метод Эйлера. Результаты вычислений приведены на рис. 2.12–2.13. Концентрация, мольн.доли 1 0,8 0,6 0,4 0,2 2 4 н‐С8Н18, время контакта 6с. н‐С8Н18, время контакта 3с. 6 8 10 Время, с и‐С8Н18, время контакта 6с. и‐С8Н18, время контакта 3с. Рис. 2.12. Зависимость концентраций реагирующих веществ от времени 70 Температура , К 620 600 580 560 2 4 6 8 10 Время, с время контакта 3с время контакта 6с Рис. 2.13. Зависимость изменения температуры от времени На основании полученных результатов можно судить об изменении концентрации веществ и температуры в реакторе идеального смешения, рассчитать степень превращения компонентов. 2.5.5. Исследование химического процесса, протекающего в реакторе идеального вытеснения в стационарном режиме Исследование закономерностей протекания химической реакции в реакторе идеального вытеснения методом математического моделирования заключается в определении концентраций реагирующих веществ на выходе из реактора и температуры потока в зависимости от времени контакта. Пусть в реакторе идеального вытеснения (РИВ) протекает химическая реакция k1 k2 A   B  C  D . Так как в реакторе идеального вытеснения состав реагентов и температура потока изменяются по длине (или времени контакта) аппарата, процесс в нем описывается системой дифференциальных уравнений (2.88, 2.89). 71 Тогда математическая модель химического процесса может быть записана в виде следующей системы уравнений материального и теплового балансов (режим работы реактора – стационарный): dC A  k1  CА ; d dCB  k1  C A  k2  CB ; d dCC  k 2  CB ; (2.95) d dCD  k 2  CB ; d dT (Q1  k  C A  Q2  k2  CB )  R '  T / p  , d Cp где k1, k2 – константы скоростей реакций; СA, СB, СC, СD – концентрации компонентов. Для решения системы дифференциальных уравнений использован численный метод Эйлера. Результаты вычислений приведены на рис. 2.14–2.15 С, мольн.доли 1 0,8 0,6 0,4 0,2 2 4 6 8 10 τ, c Рис. 2.14. Изменение концентрации компонентов в реакторе идеального вытеснения от времени контакта: и-С8Н18 н-С8Н18 72 С4Н10, С4Н8 Т, К 630 610 590 570 550 530 2 4 6 8 10 , c Рис. 2.15. Зависимость изменения температуры в реакторе идеального вытеснения от времени контакта X, % 100 80 60 40 20 2 4 6 8 10 , c Рис. 2.16. Зависимость степени превращения от времени контакта Вопросы для самоконтроля 1. 2. Какие конструкции гомогенных реакторов применяются в химической технологии? Дайте классификацию химических реакторов. 73 3. 4. Приведите примеры гомогенных химических процессов. Какие гидродинамические модели потоков наиболее широко применяются при моделировании химических реакторов? 5. В чем состоит сущность иерархического построения математической модели химического реактора? 6. каково практическое применение результатов математического моделирования химических реакторов? 7. Какими системами уравнений описываются математические модели гомогенных химических реакторов? 8. Какие численные методы применяются для исследования математических моделей гомогенных химических реакторов? 9. Назовите принципы построения математических моделей изотермических реакторов: идеального смешения, идеального вытеснения. 10. Охарактеризовать уравнения теплового балансов реакторов: адиабатический и политропический режимы работы. 74 3. ЭКСПЕРИМЕНТАЛЬНО-СТАТИСТИЧЕСКИЕ МЕТОДЫ ПОСТРОЕНИЯ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ 3.1. Основные понятия и определения В общем случае при моделировании химико-технологических процессов (ХТП) необходимо знание физико-химических закономерностей их протекания и экспериментальных данных для проверки адекватности моделей [1]. Однако далеко не всегда имеется возможность детального изучения механизма и физико-химической сущности химико-технологических процессов. В то же время, задачу оптимизации и управления такими процессами решать необходимо. В этих случаях разрабатывают так называемые эмпирические модели с применением экспериментально-статистических методов: при неизвестном механизме протекающих в объекте процессов изучают зависимость отклика системы на изменение входных параметров. В отличие от физико-химических моделей в них не учитываются закономерности протекания реальных процессов и их построение базируется на формализованном описании экспериментальных данных [1, 25]. Математическое описание объекта в этом случае будет представлять собой систему эмпирических зависимостей, полученных в результате статистического обследования объекта. Эти модели называются статистическими и имеют вид корреляционных или регрессионных соотношений между входными и выходными параметрами объекта. Естественно, в структуре уравнений статистических моделей не отражены физические свойства объекта моделирования. Основной и необходимый источник информации для построения статистической модели – эксперимент, а обработка экспериментальных данных осуществляется методами теории вероятности и математической статистики. Технологический объект в этом случае представляется в виде «черного ящика» (рис. 3.1). Математической моделью объекта будет функция отклика y    x1 , x2 , xn , b1 ,..., bn  , 75 (3.1) где y – выходной параметр процесса; x1,…,xn – независимые переменные, которые варьируются при постановке эксперимента; b1 ,..., bn – коэффициенты эмпирической модели. Рис. 3.1. Схематическое изображение объекта Конкретный вид функциональной зависимости (3.1) и значения коэффициентов определяются из опытных данных. В дальнейшем будем называть:  факторами – независимые переменные x1,…, xn;  факторным пространством – пространство с координатами x1,…, xn;  поверхностью отклика – геометрическое изображение функции отклика в факторном пространстве. В том случае, когда исследование поверхности отклика ведется при неполном знании механизма изучаемых явлений, аналитическое выражение функции отклика неизвестно, поэтому математическая модель представляется в виде полинома N N N i 1 i , j 1 i j i 1 y   0    i xi    ij xi x j    ii  xi 2   , (3.2) где βi, βij, βii – теоретические коэффициенты, характеризующие соответственно линейные эффекты, эффекты взаимодействия и квадратичные эффекты. Они называются коэффициентами регрессии, а уравнение (3.2) – уравнением регрессии. Коэффициенты регрессии:   1  ; 2  ;; x1 x2 1,2  2  ,; x1x2  2  2 . 11  ,  22  2x12 2x2 2 76 Результат эксперимента на сложном – объекте обычно величина случайная. Это может быть обусловлено погрешностью измерений, иногда случайными воздействиями («шумами»). Значения выходных измерений, как правило, отличаются друг от друга. Поэтому при обработке экспериментальных данных можно определить только так называемые выборочные коэффициенты регрессии: b0,bi,bij,bii…, которые являются лишь оценками для теоретических коэффициентов регрессии  (коэффициенты регрессии, которые можно было бы получить для некоторой генеральной совокупности, состоящей из всех мыслимых опытов). В результате пользуются приближенным уравнением регрессии, полученной по ограниченной выборке экспериментальных данных: N N N yˆ  b0   bi xi   bij xi x j   bii xi 2  , i 1 i , j 1 (3.3) i 1 где ŷ – выборочная оценка для y (предсказанное значение выходного параметра); b0 – свободный член уравнения регрессии; bi,bij,bii – коэффициенты регрессии, характеризующие соответственно линейные эффекты, эффекты взаимодействия и квадратичные эффекты. Уравнение регрессии (3.3) используется для построения статистических моделей объектов химической технологии. С точки зрения исследования физико-химических свойств процессов эта модель не несет никакой информации. Справедлива такая модель только для объекта, на котором проводили эксперимент. Однако такие модели широко используются при решении задач оптимизации. Конкретный вид эмпирических моделей (3.1) определяется по результатам экспериментов – активных или пассивных. 3.2. Статистические модели объектов на основе пассивного эксперимента Эмпирические модели строятся на основе пассивных и активных экспериментов. Эмпирическая модель представляет собой уравнение, описывающее экспериментальные данные. Известно, что измеряемые величины – случайные величины. Поэтому в соответствии с законами теории вероятности и математической статистики обработку экспериментальных данных при построении математической модели проводят с применением статистических методов [1, 5, 25]. Суть пассивного эксперимента: исследователь собирает некоторый объем экспериментальной информации, т. е. значений параметров (факторов) xi и выходного параметра yi, причем происходит это в режиме 77 нормальной эксплуатации объекта. Данные (выборка) берутся с промышленной или с лабораторной установок. В разделе 3.1 показано, что в общем виде эмпирические модели могут быть представлены в виде приближенных уравнений регрессии (3.1):   y   x1 , x2 , xn , b1 ,..., bn . (3.4) Для получения конкретного вида эмпирической модели (3.4) необходимо выполнить следующее:  найти конкретный вид функции  в уравнении (3.4);  определить значения коэффициентов регрессии bi;  выполнить статистический анализ полученных результатов. Для получения статистических математических моделей в виде полиномов используют методы корреляционного и регрессионного анализа. Процесс построения статистической модели состоит из нескольких этапов:  записывается уравнение модели в виде полинома n–й степени.  рассчитываются коэффициенты этого полинома;  оценивается наличие линейной связи между факторами, т. е. рассчитывается коэффициент парной корреляции;  оценивается значимость коэффициентов полинома по критерию Стьюдента (t);  устанавливается адекватность уравнения регрессии реальному процессу по критерию Фишера (F). 3.2.1. Методы корреляционного и регрессионного анализа Методы корреляционного и регрессионного анализов широко применяются для выявления и описания зависимостей между случайными величинами по экспериментальным данным и базируются на теории вероятности и математической статистике [1, 25]. Корреляционный анализ основывается на предпосылке о том, что переменные величины y (выходной параметр) и xi (факторы) являются случайными величинами и между ними может существовать так называемая корреляционная связь, при которой с изменением одной величины изменяется распределение другой. Для количественной оценки тесноты связи служит выборочный коэффициент корреляции. Можно выделить три типа коэффициента корреляции: 1. Простой коэффициент корреляции, или коэффициент парной корреляции, определяет величину (тесноту) зависимости между двумя переменными (x или y) и определяется по формуле 78    n  x x y y i i r  i 1 , xy  N  1  S x  S y (3.5) где x , y – среднеарифметические значения переменных x, y; N – число опытов; Sx, Sy – среднеквадратические отклонения случайных величин: Sx  N 2  xi  x i 1 ; Sy  N 1   N   y y i 1 i N 1  2 . (3.6) Коэффициент корреляции характеризует степень тесноты линейной зависимости. Если случайные величины x и y связаны точной линейной функциональной зависимостью ŷ  b0  b1 x , то r  1 , при этом знак xy коэффициента корреляции соответствует знаку коэффициента b1 . В том случае, когда величины x и y связаны произвольной стохастической зависимостью, коэффициент корреляции может принимать значение в интервале 1  r  1. (3.7) xy Если r xy  0 , x и y – независимы, корреляции нет. При r xy  0 суще- ствует положительная корреляционная связь между x и y (с ростом x увеличивается y), если r  0 – отрицательная (с ростом x уменьшается y). xy Коэффициент корреляции, существенно отличающийся от 1, характеризует слишком большую долю случайности и слишком большую долю криволинейности связи между случайными величинами. О наличии или отсутствии корреляции между двумя случайными величинами качественно можно судить по виду поля корреляции (рис. 3.2). а б Рис. 3.2. Поля корреляции случайной величины: а – сильная положительная корреляция между x и y; б – слабая корреляция; в – корреляции нет 79 в Оценка зависимости случайных величин по выборочному коэффициенту корреляции называется корреляционным анализом. При вычислении коэффициента корреляции удобно пользоваться следующими формулами: x y (3.8)  ( xi  x )( yi  y )   xi yi   Ni  i ; ( xi ) 2 2  xi  N 1 2 2 2 2 ( N  1) S x   xi  ( xi ) ; S x  ; N N 1 ( yi ) 2 2  yi  N 1 2 2 2 2 ( N  1) S y   yi  ( yi ) ; S y  , N N 1 где N – число опытов; S x2 S y2 – выборочные дисперсии величин x. 2. Коэффициент частной корреляции измеряет линейную зависимость между двумя переменными после устранения части зависимости, обусловленный зависимостью этих переменных с другими переменными. При исследовании зависимости y от x1 и x2 наличие корреляции между x1 и x2 и между y и x2 будет влиять на корреляцию между y и x1. Для того чтобы устранить влияние x2, необходимо измерить корреляцию между y и x1 при x2 = const. Частный коэффициент ryx x оценивает степень влияния фактора 1 2 x1 на y при условии, что влияние x2 на y исключено: ryx1  ryx2  rx2 x1 (3.9) ; ryx1x2  1 1 1  ryx2 2 2 1  rx21x2 2  ryx2 x1    ryx2  ryx1  rx2 x1 1 2 1  r  1  r  2 yx1 2 x1 x2 1 2 . 3. Множественный коэффициент корреляции определяет величину зависимости одной переменной от нескольких. Коэффициент корреляции показывает, существует ли связь между x и y, но самого вида функции не дает. Для характеристики формы связи при изучении корреляционной зависимости пользуются уравнением приближенной регрессии. Регрессионный анализ предполагает (рассматривает) связь между зависимой величиной y и независимыми переменными x1,…,xi. Эта связь представляется с помощью математической модели, т. е. уравнения, которое связывает зависимую и независимые переменные [1, 25]. 80 Перечислим предпосылки, на которых базируется регрессионный анализ: 1. Результаты наблюдений y1, y2,…,yn представляют собой независимые, нормально распределенные случайные величины. 2. Входные факторы xi измеряются с пренебрежимо малой ошибкой по сравнению с ошибкой измерения y. 3. Выборочные дисперсии S12 , S 22 ,..., S n2 значения выходного параметра у, полученные при одинаковом количестве параллельных опытов, должны быть однородны. Обработка экспериментальных данных при использовании корреляционного и регрессионного анализа дает нам возможность построить статистическую математическую модель в виде уравнения регрессии. Таким образом, методы корреляционного и регрессионного анализов тесно связаны между собой. Ранее отмечалось, что коэффициент корреляции дает нам меру зависимости между двумя величинами, но не дает вида (формы) взаимосвязи. Для характеристики формы связи пользуются уравнением регрессии. Постановка задачи По данной выборке объема n найти уравнение приближенной регрессии и оценить допускаемую при этом ошибку. Эта задача решается методом корреляционного и регрессионного анализов. Нужно найти yˆ  f ( x ) . По сгущениям точек можно найти определенную зависимость, т. е. получить вид уравнения регрессии. Если разброс точек значительный, то регрессии не будет. Вид уравнения регрессии зависит от выбираемого метода приближения. Обычно пользуются методом наименьших квадратов:  N F   y  f (x ) i i i 1  2  min, (3.10) или  N F   y  yˆ i i i 1  2  min, (3.11) где yi , yˆi – экспериментальные и расчетные значения выходной величины соответственно. 81 При обработке результатов пассивных экспериментов получают линейные и нелинейные эмпирические модели, которые должны достаточно точно описывать всю совокупность опытных данных. Сложность в этом случае заключается в правильном выборе вида модели и определении параметров модели (коэффициентов уравнения). Рассмотрим различные случаи приближенной регрессии [1, 5, 25]. 3.2.1.1. Линейная регрессионная модель с одной независимой переменной Рассмотрим случай линейной регрессии от одного параметра. При моделировании химико-технологических процессов (ХТП) во многих случаях связь между входными (x) и выходными (y) параметрами можно аппроксимировать линейным уравнением y  b  b x. 1 (3.12) Для получения вида математической модели необходимо определить коэффициенты уравнения регрессии b0 и b1. Для этого применим метод наименьших квадратов (МНК):  N F   y b b x i 0 1i i 1  2  min. Таким образом, процедура нахождения коэффициентов регрессии сводится к задаче определения минимума функции. Необходимым условием экстремума функции является равенство нулю частных производных функции по искомым величинам (коэффициентам): N  F 2     yi  b0  b1xi   1  0;  b i 1  F N   2   yi  b0  b1 xi   xi  0; i 1  b1    yi  b0  b1 xi   0;     yi  b0  b1 xi   xi  0; (3.13)  nb0  b1  xi   yi ;  2 b0  xi  b1  xi   xi yi . Решаем систему уравнений методом наименьших квадратов. В результате получаем формулы для вычисления коэффициентов b0 и b1: 82 b0   yi  xi  xi yi  xi N  xi  xi  xi 2 2   yi  xi   xi yi  xi N  xi    xi  2 2 2 ;  yi N  xi yi   xi  yi  xi  xi yi b1  .  2 N  xi N  xi2    xi  2  xi  xi (3.14) N (3.15) После вычисления коэффициентов уравнения (3.12) приступают к исследованию полученной математической модели. Такое исследование называется статистическим (регрессионным) анализом. Рассмотрим последовательность данного анализа. 3.2.1.2. Статистический анализ результатов 1. Для оценки тесноты линейной зависимости между факторами рассчитывают коэффициенты парной корреляции rxy по формуле (3.5). Чем ближе значение rxy к 1, тем вероятнее наличие линейной связи. Следовательно, зависимость между x и y в определенном диапазоне будет иметь вид yˆ  b0  b1 x. 2. Проверка однородности дисперсий. 1) Определяется среднее по результатам параллельных опытов (если есть параллельные опыты): m  yiu yi  u 1 ; m i  1,..., N , (3.16) где m – число параллельных опытов; N – количество опытов в выборке. 2) Определяются выборочные дисперсии: m Si2    yiu  yi  u 1 m 1 2 ; 3) Суммируются дисперсии: N 2  Si ; i 1 83 i  1, N (3.17) 4) Выбирается максимальная дисперсия, составляется отношение 2 Smax G N S i 1 (3.18) , 2 i 2 где S max – максимальное значение выборочной дисперсии. Проверяется однородность дисперсий по критерию Кохрена (G) (при одинаковом количестве параллельных опытов). (q, f), то дисперсии однородны, где q – уровень знаЕсли G  G табл чимости; f – число степеней свободы. Число степеней свободы f1 = m – 1; f2 = N. 5) Определяется дисперсия воспроизводимости N 2 i 1 S воспр  3. 2  Si N ( m  1) . (3.19) Число степеней свободы f = N(m–1) – для одинакового числа опытов m. Оценивается значимость коэффициентов полинома по критерию Стьюдента (t) (предпосылка – отсутствие корреляции между факторами): tbi  bi , Sbi (3.20) где bi – i-й коэффициент уравнения регрессии; Sbi – среднеквадратическое отклонение i-го коэффициента. Для случая линейного полинома Sb0 и Sbi вычисляются по следующим формулам: N Sb0  2 Sвоспр  xi i 1   N N N  xi2   i i 1 S  b1 2 ; (3.21) . (3.22) i 1 2 Sвоспр N N   N N  x   xi i 1 2 i 2 i 1 При числе факторов больше двух подобные формулы из-за громоздкости не используются. Если t  t ( q , f ) , то коэффициент bi значим (значимо отличаетb i табл ся от 0). В противном случае – не значим. 84 4. Проверка модели (полученного уравнения) на адекватность осуществляется по критерию Фишера (F). Полагают, что уравнение регрессии адекватно описывает исследуемый процесс, если остаточная дисперсия выходной величины y, рассчитанной по уравнению регрессии относительно экспериментальных данных, не превосходит ошибки опыта. Если 2 Sост F  2  FT  q, f1 , f 2  , (3.23) Sвоспр то модель адекватна (т. е. линейное уравнение регрессии адекватно описывает исследуемый объект). Для одинакового числа параллельных опытов m1  m2  ...mn выражение для остаточной дисперсии имеет вид N  2   yi  yˆi  2  i 1 Sост , (3.24) N l где yi – среднее значение выходного параметра по результатам параллельных опытов (3.16); yˆi – расчетное значение выходного параметра. Если при проведении эксперимента опыты не дублировались, то выражение будет следующим: N   yi  yˆi  2 2  i 1 Sост , (3.25) N l где f1  N  l и f 2  N  m  1 – число степеней свободы числителя и знаменателя соответственно; l  n  1 – число членов аппроксимирующего полинома (число коэффициентов регрессии, включая свободный член); N – общее количество опытов; n – количество факторов (x1,x2,…); y – экспериментальное значение выходного параметра. i Если при проведении эксперимента не было возможности выполнить параллельные опыты, то вместо проверки модели на адекватность выполняется оценка качества аппроксимации уравнением. Это достига2 ется сравнением остаточной дисперсии Sост с дисперсией относительно среднего S y2 : N S y2    yi  y  i 1 2 , N 1 где yi – экспериментальное значение выходного параметра. 85 (3.26) 1 N  yi – среднее значение выходного параметра. N i 1 По критерию Фишера S y2 F 2 . (3.27) Sост В этом случае чем больше значение F превышает табличное Fтабл(q,f1,f2), тем уравнение регрессии эффективнее для выбранного уровня значимости q и чисел степеней свободы f1  N  l и f 2  N  l . Таким образом, критерий Фишера показывает, во сколько раз уменьшается рассеивание относительно полученного уравнения регрессии по сравнению с рассеиванием относительно среднего. y 3.2.1.3. Параболическая регрессионная модель При составлении статистических моделей объектов химической технологии часто возникает необходимость использовать уравнения нелинейной формы: чаще всего – параболу второй, третьей и более высоких степеней. В таких случаях используют метод регрессионного анализа для составления статистической модели в виде полинома второй (или более высокой) степени: N y  b   b x   b x x   b x 2  ...,... i i ij i j ij i i 1 i  j Дано уравнение y  b0  b1x  b2 x 2 , требуется определить b0, b1, b2. Коэффициенты регрессии определяют по МНК:  F   y  y  2  min ; (3.28) N F    yi  b0  b1  xi  b2  xi2   min ; i 1 2       N 2  2   yi  b0  b1  xi  b2  xi  1  0; i 1 b0 N F 2  2   yi  b0  b1  xi  b2  xi  xi  0; i 1 b1 N F 2 2  2   yi  b0  b1  xi  b2  xi  xi  0; i 1 b2 F 86 N N 2 N  b0  N  b1   xi  b2   xi   yi ;  i 1 i 1 i 1  N N 2 N 3 N  b0   xi  b1   xi  b2   xi    xi  yi ; i 1 i 1 i 1 i 1  N 2 N 3 N 4 N 2 b   x  b   x  b   x   x  y . i 2 i 1 i i 1 i  0 i 1 i 1 i 1 i Введем обозначения:   N N N 2 N 4 3 S1   xi ; S 2   xi ; S3   xi ; S 4   xi ; i 1 i 1 i 1 i 1 N N N 2 S5   yi ; S6   xi  yi ; S7   xi  yi . i 1 i 1 i 1     С учетом принятых обозначений система будет иметь следующий вид:  b0  N  b1  S1  b2  S2  S5 ;   b0  S1  b1  S2  b2  S3  S6 ; b  S  b  S  b  S  S . 0 2 1 3 2 4 7 Неизвестные коэффициенты b0, b1, b2 определяются по следующим выражениям: S5 S1 S 2 S 6 S 2 S3 S S3 S 4 S5 S 2 S 4  S6 S3 S2  S7 S1S3  S7 S2 S 2  S6 S1S4  S5 S3 S3 b0  7 ;  N S1 S 2 NS 2 S 4  S1S3 S2  S2 S1S3  S2 S2 S 2  S1S1S4  NS3 S3 S1 S 2 S3 S 2 S3 S 4 n S5 S 2 S1 S6 S3 S S7 S4 NS6 S 4  S1S7 S2  S2 S5 S3  S2 S6 S2  S1S5 S4  NS7 S3 (3.29) b1  2  ; N S1 S2 NS 2 S 4  S1S3 S2  S 2 S1S3  S2 S2 S2  S1S1S4  NS3 S3 S1 S2 S3 S 2 S3 S 4 N S1 S5 S1 S 2 S6 S S3 S7 NS2 S7  S1S3 S5  S2 S1S6  S2 S2 S5  S1S1S7  NS3 S6 b2  2 .  N S1 S 2 NS2 S4  S1S3 S 2  S2 S1S3  S2 S2 S2  S1S1S4  NS3 S3 S1 S 2 S3 S 2 S3 S 4 87 Решая систему уравнений, вычисляем коэффициенты b0, b1, b2. Аналогично будут определяться коэффициенты параболы любого порядка. Исследование уравнения проводится по статистическим критериям, так же как и в случае линейной регрессии. Адекватности уравнения регрессии эксперименту можно добиться, повышая степень полинома. Однако при этом все коэффициенты следует вычислять заново (т. к. существует корреляция между коэффициентами). Таким образом, пассивные методы сбора экспериментальной информации имеют определенные преимущества, которые заключаются в том, что информация собирается в режиме нормальной эксплуатации объекта. Это преимущество способствовало широкому внедрению регрессионного анализа для изучения объектов химической технологии. Однако полученные на базе пассивного эксперимента модели во многих случаях оказываются неэффективными. Причиной является не сам метод, а то, что не выполняются основные предпосылки регрессионного анализа: факторы измеряются с большими ошибками. В пассивном эксперименте, как правило, ошибка измерения x соизмерима, а то и больше ошибки при измерении y. Иногда ошибка измерения входных параметров превышает даже интервал изменения самих факторов. Кроме того, факторы (xi) или коэффициенты bi имеют между собой корреляционную связь. Это затрудняет статистический анализ и интерпретацию результатов. Вопросы для самоконтроля 1. В каких случаях прибегают к построению статистических моделей? 2. На чем базируется построение статистических моделей? 3. Каков общий вид статистических моделей? 4. Приведите два вида эксперимента, используемые для построения статистических моделей. 5. В чем разница между пассивным и активным экспериментом? 6. Для чего проводят корреляционный анализ? 7. Какова основная характеристика корреляционного анализа? 8. Какова суть регрессионного анализа? 9. Перечислите виды регрессии, привести примеры. 10. Назовите метод, применяемый для оценки коэффициентов уравнения регрессии. 11. Какова последовательность статистического (регрессионного) анализа? 88 3.3. Статистические модели на основе активного эксперимента (методы планирования экстремальных экспериментов) Обработка результатов пассивного эксперимента проводится методами корреляционного и регрессионного анализов, и выбор вида эмпирической модели является достаточно сложной задачей, т. к. вид уравнения регрессии необходимо определить по характеру изменения переменных на графике эмпирической линии регрессии, полученной по выборке экспериментальных данных. Стремление избавиться от недостатков моделей на базе пассивного эксперимента привело к математической теории активного эксперимента [26]. Активный эксперимент ставится по заранее составленному плану и обрабатывается по некоторому оптимальному алгоритму с целью составления математической модели. Одним из основных методов теории активного эксперимента является статистическое планирование эксперимента [1, 5, 25, 26]. 3.3.1. Планы первого порядка 3.3.1.1. Полный факторный эксперимент При планировании по схеме полного факторного эксперимента (ПФЭ) реализуются все возможные комбинации факторов на всех выбранных для исследования уровнях. Суть факторного эксперимента:  одновременное варьирование всех факторов при проведении эксперимента по определенному плану;  представление математической модели (функции отклика) в виде линейного полинома;  исследование полученного полинома методами математической статистики. Необходимое количество опытов N при ПФЭ определяется по формуле N = ln, (3.30) где n – число факторов; l – число уровней, на которых варьируются факторы. Уровни факторов – это границы исследуемой области по данному технологическому параметру. Выбор экспериментальной области по каждому фактору связан с тщательным анализом экспериментальной информации. Обычно применяется планирование на двух уровнях, т. е. l = 2, тогда при n = 2 число опытов будет равно N = 22 = 4. 89 Основной (нулевой) уровень (центр плана эксперимента) – это некоторое начальное значение фактора при составлении математической модели, это точка с координатами ( x10 ,..., xn0 ) . Интервал варьирования – часть области определения фактора, симметричная относительно его основного уровня ( x ). Пример. На процесс влияют два фактора: температура (x1) и концентрация вещества (x2). Известно априори, что температура изменяется в интервале 100–200 °С; концентрация 20–30 %. Для двух факторов, без учета их взаимодействия, соответствующая модель может быть записана как ŷ  b0  b1 x1  b2 x2 . Приведем условия эксперимента: Верхний уровень: Нижний уровень: Основной уровень: Интервал варьирования Основной уровень рассчитывается как x10  x1, °С 200 100 150 50 x2, % 30 20 25 5 x1min  x1max 200  100   150; 2 5 x2min  x2max 30  20 x    25. 2 2 Интервалы варьирования: 2 x1max  x1min 200  100 x1    50; 2 2 x2max  x2min 30  20 x2    5. 2 2 В координатах на плоскости это можно представить (рис. 3.3): План эксперимента указывает расположение в n-мерном пространстве опытных точек независимых переменных или условия всех опытов, которые необходимо провести. При ПФЭ эксперимент ставится только на границе области, т. А – центр области. В большинстве случаев эксперимент задается в виде матрицы планирования – это план (таблица), каждая строчка которого 90 представляет собой условия опыта, а каждый столбец матрицы соответствует значениям переменных в различных опытах. Рис. 3.3. Схема расположения опытных точек полного факторного эксперимента 22 Составим матрицу планирования для предыдущего примера. Количество опытов в матрице планирования N = 22 = 4 (табл. 3.1). Таблица 3.1 Матрица планирования N 1 2 3 4 x2 30 20 30 20 x1 200 200 100 100 Y Y1 Y2 Y3 Y4 Это матрица планирования в натуральном масштабе. Матрица планирования составляется для того, чтобы эксперимент провести по определенному плану, определить значения выходного параметра в каждом опыте и на основании экспериментальных данных построить статистическую модель. При планировании первого порядка получают математическую модель вида yˆ  b0  b1 x1  b2 x2  ...  bn xn . Кодирование переменных Для удобства расчетов перейдем от натуральных координат (натуральных единиц измерения) к безразмерным. Формула перехода (или кодирования) имеет вид xi  xi0 Xi  , (3.31) xi 91 где xi – значения (верхний или нижний уровень) натуральной переменной; xi0 – основной уровень натуральной переменной; xi – интервал варьирования натуральной переменной; Xi – кодированное значение i-го фактора (на верхнем или на нижнем уровне). Перейдем от натуральных переменных к кодированным: для температуры: 200  150  1; 50 100  15 н x1  50  1. x в 1  для давления: 30  25  1; 5 20  25 н x 2  5  1. x в 2  Фактически мы обозначили значения факторов на верхнем уровне +1, (200, 30), а на нижнем (100, 20) –1; Представим матрицу планирования в кодированных единицах. N x0 x1 x2 1 2 3 4 +1 +1 +1 +1 +1 +1 –1 –1 +1 –1 +1 –1 N x0 x1 x2 1 2 3 4 + + + + + + – – + – + – или Это матрица планирования в безразмерном масштабе. x0 – фиктивная переменная (+1), необходимая для вычисления свободного члена полинома. 92 Расположение опытных точек в факторном пространстве будет следующим: Свойства матрицы планирования Матрица планирования в безразмерных единицах обладает следующими оптимальными свойствами:  ортогональностью: скалярное произведение двух любых столбцов матрицы равно нулю: N  xui x ji  0; uj; u, i = 1,…,n; i 1 (3.32)  симметричностью: сумма элементов всех столбцов матрицы, кроме первого, равна нулю: N  x  0, u = 1,…,n; iu i 1 (3.33)  нормировкой: сумма квадратов элементов каждого столбца равна числу опытов N 2  xiu  N , u = 1,…,n; i 1 (3.34)  свойство ротатабельности: все точки в матрице планирования подбираются так, что точность предсказания значений выходного параметра одинакова на равных расстояниях от центра эксперимента и не зависит от направления. На основании всех перечисленных выше свойств, в частности ортогональности и ротатабельности, значительно упрощается расчет коэффициентов регрессии. 93 Расчет коэффициентов регрессии После того как составлен план (матрица планирования), проводят эксперименты (дублируя опыты) и на основании результатов рассчитывают коэффициенты в уравнении регрессии по формулам: 1 N  x0i  yi ; N i 1 1 N bi   xiu  yi ; (3.35) N i 1 1 N bij   xiu xij yi . N i 1 Приведенные формулы получены на основании метода наименьших квадратов (см. 3.14, 3.15) с учетом свойств матрицы планирования (3.32–3.34). Пример. Вычислить коэффициенты регрессии на основании экспериментальных данных, приведенных в табл. 3.2. b0  Таблица 3.2 N x0 x1 x2 1 2 3 4 +1 +1 +1 +1 1 1 –1 –1 1 –1 1 –1 x1x2 Y 1 –1 –1 1 68,2 62,5 56,8 50,7 68, 2  62,5  56,8  50,7  59,6; 4 68, 2  62,5  56,8  50,7 b1   5,8; 4 68, 2  62,5  56,8  50,7 b2   2,95; 4 68, 2  62,5  56,8  50,7  0,1. b12  4 С учетом полученных коэффициентов запишем уравнение регрессии: yˆ  59,6  5,8 x1  2,95 x2  0,1x1 x2 . b0  После вычисления коэффициентов регрессии приступают к статистическому анализу уравнения регрессии. 94 3.3.1.2. Статистический анализ уравнения регрессии 1. Статистический (регрессионный) анализ состоит из следующих этапов: Оценка дисперсии воспроизводимости: 1) Определяются средние значения выходного параметра по результатам параллельных опытов: N = ln; (3.36) m yi   yiu u 1 (3.37) i  1,..., N , ; m где m – число параллельных опытов. 2) Рассчитываются выборочные (построчные) дисперсии: m 2 Si    yiu  yi  2 u 1 m 1 (3.38) . N 3) Суммируются построчные дисперсии  S i 2 , находится максиi 1 2 мальное значение дисперсии S max . 4) Находится отношение G S max S 2 . 2 (3.39) i Проверяется дисперсия на однородность по критерию Кохрена (G), если G5. Если линейное уравнение регрессии в результате статистического анализа оказалось неадекватным, то поступают следующим образом:  добавляют к ядру плана (ПФЭ или ДФЭ) 2n точек, расположенных на осях координат факторного пространства на расстоянии   от центра плана. Эти точки называются «звездными». Координаты «звездных» точек ( , , . Величина  – «звездное плечо».  увеличивают число экспериментов в центре плана n0: точки с координатами (0, …, 0). На рис. 3.4 приведена схема ЦКП для n = 2. Точки 1234 – ПФЭ 22; точки 5678 – «звездные» точки с координатами (   и (0,  ). Количество опытов в матрице ЦКП определяется как N = 2n + 2n + n0, если n5; (3.45) N = 2n–1 + 2n + n0, если n>5. Построим композиционный план для n = 2 (табл. 3.4). Количество опытов N = 22 + 22 + 1 = 9.      Рис. 3.4. Схема ЦКП для n = 2 100 (3.46) Таблица 3.4 x1 + – + – + – x0 + + + + + + + + + N 1 2 3 4 5 6 7 8 9 x12 + + + + x2 + + – – + – 2 2 x22 + + + + 2 2 Эта матрица не ортогональная, т. к. x x 0j x 2ij  0; x  0. 2 ij uj Ортогональность композиционных планов достигается выбором значения «звездного плеча» . Приведем некоторые значения  для ортогональных планов (n0 = 1). 2(22) 1.0 n  3(23) 1.215 4(24) 1.414 5(25) 1.547 Уравнение регрессии при ортогональном ЦКП в общем виде будет следующим (например, для двух факторов): ŷ  b0  b1 x1  b2 x2  b12 x1 x2  b11 x1  b22 x2 . (3.47) Величины x1* и x2* введены для того, чтобы привести матрицу планирования к ортогональному виду, а коэффициенты bi определялись независимо друг от друга:  x ji  x ji 2 1 N 2   x ji , N j 1 (3.48) где j – номер опыта, i – номер фактора. Для того чтобы получить уравнение регрессии в обычной форме: ŷ  b0  b1 x1  b2 x2  b12 x1 x2  b11 x12  b22 x2 2 , (3.49) находят величину b0  b0  b11 N b x j12  22  x j 2 2 .  N j 1 N 101 (3.50) Приведем матрицу ортогонального ЦКП для n = 2,  = 1 (табл. 3.5). Таблица 3.5 N x1 x2 x1x2 x1 x2 1 2 3 4 5 6 7 8 9 +1 –1 +1 –1 +1 –1 +1 +1 –1 –1 +1 –1 +1 –1 –1 +1 0,33 0,33 0,33 0,33 0,33 0,33 –0,67 –0,67 –0,67 0,33 0,33 0,33 0,33 –0,67 –0,67 0,33 0,33 –0,67 y Значения x1 и x2 рассчитываются по формуле (3.47). Например: x11*  (1) 2  6  1  0,67  0,33; 9 6  0,33; 9 *  0  0,67  0,67. x52 * 1 x51 Для пересчета значений факторов в натуральные единицы пользуемся формулой пересчета xi  xi 0 Xi  ; xi . xi  X i  xi  xi 0 . Матрица (3.5) ортогональная, т. е. x x 2  0;  xij 2 xuj 2  0 , 0 j ij но не ротатабельная. Здесь j – номер опыта. Коэффициенты регрессии при ортогональном ЦКП рассчитываются по следующим формулам: 102 1 N b0   y j ; N j 1  N x j 1 N bi  ij yj  x  ; 2 ij j 1 N biu  x x j 1 N ij uj N bij  ;  (x x j 1 x  ij j 1 N (3.51) yj ij uj ) 2 yj  x  2 .  ij i 1 Регрессионный анализ уравнения проводится по схеме, приведенной ранее. Для расчета дисперсий при определении коэффициентов регрессии используют выражения S 2 bo Sвоспр 2  N Sb0 2  S 2b   ; nSbij N Sbi  Sвоспр 2 2 N  x  j 1 Sbiu   x N ji x ju   x  j 1 bj Sbj ; ; S 2 воспр 2 2 ji ji N j 1 ti  2 Sвоспр 2 2 Sbii  x 2  ji , 103 ; 2 ; (3.52) где i≠u, i≠0. Коэффициент значим, если tj>tT(q,f2); f2 – число степеней свободы 2 S воспр. Заключительный этап – проверка уравнения на адекватность по критерию Фишера. Ротатабельные планы второго порядка [25] Ротатабельные планы были предложены в 1957 г. Боксом и Хантером. Этот метод планирования эксперимента позволяет получить более точное математическое описание поверхности отклика по сравнению с ортогональным ЦКП. Это достигается за счет увеличения числа опытов в центре плана и специального выбора величины звездного плеча . Приведем некоторые значения  и n0 для различного числа факторов n. Параметр плана Количество факторов 2 3 4 5 5 6 6 7 7 22 23 24 25 25–1 26 26–1 27 27–1 1,414 1,682 2,0 2,380 2,0 2,830 2,380 3,360 2,830 5 6 7 10 6 15 9 21 14 Ядро плана  n0 Составим матрицу ротатабельного планирования второго порядка для n = 2,  = 1,414, n0 = 5 (табл. 3.6). Таблица 3.6 Nn x0 x1 x2 х1х2 x12 х22 1 2 3 4 5 6 7 8 9 10 11 12 13 + + + + + + + + + + + + + + + +1,414 –1,414 + + +1,414 –1,414 + + + + + + 2 2 + + + + 2 2 104 Матрица ротатабельного планирования второго порядка неортогональна, т. к. N x j 1 N x j 1 x 2  0; 0 j ij ij 2 xuj 2  0. Формулы для расчета коэффициентов имеют вид n 2 AB     b0  S B n 2 C S    ii ; N  i 1  C  Si ; N C 2 Sij b ji  ; i  j; BN n AC   bii   SiiC  B  n  2   n   C 1  B   Sii  2 BS0  ; N  i 1  А, В, С – константы, которые определяются как 1 A ; 2 B  n  2  B  n  bi  B n N ;  n  2  N  N0  C N , N  N0 n – число факторов; N – общее число опытов ротатабельного ЦКП; N0 –число опытов в центре плана. По результатам эксперимента вычисляют суммы: N S0   y j ; j1 N Si   x ji y j ; N i  1,..., n; Siu   xij x jk y j ; i  u; j 1 N Sii   x 2 ji y j . 105 (3.53) Оценки дисперсий в определении коэффициентов вычисляются по следующим формулам: 2 AB  n  2  2 S воспр ; N S 2b0  S 2 iu  S S 2 bii  S 2воспр N  N0 2 biu  AC 2 S 2 воспр N ;  i  1,..., n  ; C 2 S 2 воспр N ;  B  n  1   n  1  . Коэффициенты значимы, если bi  Sbi t : N S 2 ост   y э j  y j p   S 2 воспр  N 0  1  n  2  n  1  N 2 f ост  N   n  2  n  1  2  N 0  1 ; (3.53)  N 0  1. Проверку адекватности уравнения регрессии проводят с помощью критерия Фишера. После того, как получено уравнение регрессии второго порядка, адекватно описывающее почти стационарную область, его исследуют для выбора оптимальных условий технологического процесса. Полученное уравнение дает информацию о форме поверхности отклика. Для изучения конфигурации поверхности отклика уравнение приводят к канонической форме (эллиптический параболоид, седло, и т. д.) и исследуют на локальный экстремум. Вопросы для самоконтроля 1. 2. 3. 4. 5. В каких случаях приступают к планированию второго порядка? Назовите виды планов второго порядка. Какова суть ортогонального планирования второго порядка? Суть ротатабельных планов второго порядка. Ваши действия после получения полинома второго порядка? 106 3.4. Симплексный метод планирования и оптимизации Симплексный метод используется на стадии восхождения по поверхности отклика [1]. Симплексом называется правильный многогранник, имеющий n + 1 вершину, где n – число факторов, влияющих на процесс. Так, если n = 1, то симплексом является отрезок прямой, при n = 2 – правильный треугольник, при n = 3 – тетраэдр и т. д. Метод последовательного симплекс–планирования состоит в следующем: начиная восхождения, планируют исходную серию опытов так, чтобы точки, соответствующие условиям проведения этих опытов, образовывали правильный симплекс в многомерном факторном пространстве. Начальная серия опытов соответствует вершинам исходного симплекса (рис. 3.5, точки 1 2 3). Х1 Рис. 3.5. Схема движения к оптимуму Условия первых опытов выбираются из тех значений факторов, которые соответствуют наиболее благоприятным из известных технологических режимов процесса. Проводят первую серию опытов (рис. 3.5, точки 1, 2, 3), после чего выявляют точку (опыт), которая дала наихудший результат (сравнивая точки 1, 2 и 3 рис. 3.5 значения выходного параметра). На нашем рисунке это точка 1. Эту «плохую» точку заменяют новой (рис. 3.5, т. 4), представляющей собой зеркальное отображение относительно противоположной грани симплекса (т. 2 – т. 3 рис. 3.5). В новой точке (опыт 4) проводят эксперимент. Далее сравнивают между собой результаты опытов в вершинах нового симплекса (2,3,4), отбрасывают самый «неудачный» и переносят эту вершину симплекса (т. 3) в т. 5 рис. 3.5. Получают новый симплекс (2, 4, 5 рис. 3.5) и т. д. 107 Эта процедура повторяется до тех пор, пока не будет достигнут оптимум. Вершины (условия опытов) исходного симплекса задаются при помощи специальной таблицы (табл. 3.7). Условия каждого нового опыта в отраженной точке рассчитываются по формуле xi n  2  2 xi c  xi  , i  1,..., n, (3.54) где xi – значение i-го фактора в наихудшей точке (опыте) предыдущего симплекса; xic – координаты центра противоположной грани, которые находятся как n 1 xi  c x j 1 ij . (3.55) n Построение матрицы исходного симплекса Прежде чем начать движение по поверхности отклика, необходимо определить условия опытов в исходном симплексе. Для вычисления этих значений пользуются матрицей опытов исходного симплекса в кодированных переменных. Таблица 3.7 N 1 2 3 4 5 6 7 x1 0,5 –0,5 x2 0,289 0,289 –0,578 x3 0,204 0,204 0,204 –0,612 x4 0,158 0,158 0,158 0,158 –0,632 x5 0,129 0,129 0,129 0,129 0,129 –0,645 x6 0,109 0,109 0,109 0,109 0,109 0,109 –0,654 Приступая к оптимизации, необходимо при помощи таблицы рассчитать матрицу исходной серии опытов в натуральных единицах по следующим формулам: xi  x0i ; Xi  xi xi  x 0i  xi X , где Xi – кодированные значения из таблицы. 108 (3.56) При шаговом восхождении по поверхности возможны следующие случаи: 1. В результате отображения некоторой наихудшей точки симплекса в новом симплексе отраженная точка тоже оказалась наихудшей. В этом случае следует вернуться в предыдущий симплекс и двигаться из него, отбросив точку, показавшую второе наихудшее значение y. 2. Симплекс вращается вокруг некоторой точки, отвечающей наибольшему значению y. После проведения n + 1 опыта необходимо прекратить движение и повторить точку (опыт), вокруг которой вращались. Если значение в этой точке подтверждается, то, следовательно, достигнута область оптимума. Следует отметить, что симплексный метод – локальный метод поиска экстремума. При использовании симплекс-метода дублировать опыты не обязательно, т. к. ошибка в отдельном опыте может только несколько замедлить оптимизацию. 3.4.1. Пример поиска оптимальных условий методом симплекс–планирования Исследовали процесс механического обезвоживания торфа. Ставится задача: получить торф влажностью W = 60 . Факторами, влияющими на удаление влаги из торфа, являются: кг x1(g) – удельная нагрузка фильтра торфом, 2 ; м x2() – продолжительность отжатия, с; x3(p) – давление прессования, а; x4(T) – температура, °С. Сформируем условия опытов и шаги варьирования (n = 4). x3(p) x4(T) 0,3 x2() 60 1,2 60 0,2 30 0,8 30 0,5 0,1 90 30 2,0 0,4 90 30 x1(g) xi0 xi Верхний уровень Нижний уровень Количество факторов n = 4, следовательно, количество опытов в исходном симплексе n + 1 = 5. Для расчета условий опытов в исходном симплексе используем формулу кодирования (3.56) и матрицу исходного симплекса в кодах. 109 I. Значение первого фактора в пяти опытах: x11 = 0,3 + 0,20,5 = 0,4; x12 = 0,3 + 0,2(–0,5) = 0,2; x13 = 0,3 + 0,2(0) = 0,3; x14 = 0,3 + 0,20 = 0,3; x15 = 0,3 + 0,20 = 0,3. II. Значение второго фактора в пяти опытах: x21 = 60 + 300,289 = 68,7; x22 = 60 + 300,289 = 68,7; x23 = 60 – 300,578 = 42,7; x24 = 60 + 300 = 60,0; x25 = 60,0. III. Значение третьего фактора в пяти опытах: x31 = 1,2 + 0,80,204 = 1,36; x32 = 1,36; x33 = 1,36; x34 = 1,2 – 0,80,612 = 0,71; x35 = 1,2. IV. Значение четвертого фактора в 5 опытах: x41 = 60 + 300,158 = 64,71; x42 = 64,7; x43 = 64,7; x44 = 64,7; x45 = 41,0. Заполним таблицу (табл. 3.8). После расчета условий опытов в исходном симплексе реализуют пять опытов (4 + 1). Выбирают «наихудшую точку» (табл. 3.8, т. 3) и находят ее зеркальное отображение. Рассчитывают координаты отображенной точки по формулам (3.54) и (3.55). Для этого суммируют значения xi, кроме значений в т. 3 (наихудшая): 0, 4  0, 2  0,3  0,3 x1c   0,3; 4 2  68,7  60  2 x2 c   64,39; 4 2  1,36  0,72  1, 2 x3c   1,16; 4 64  7  3  41 x4 c   58,8. 4 110 Таблица 3.8 N x1 x2 x3 x4 y5(W,) 1 2 3 4 5 6 7 8 9 10 11 12 0,4 0,2 0,3 0,3 0,3 0,3 0,3 0,3 0,3 0,15 0,176 0,12 68,7 68,7 42,6 60,0 60,0 86,2 81,8 92,7 76,3 93,3 94,1 88,2 1,36 1,36 1,36 0,72 1,2 0,96 1,72 1,5 0,87 0,98 1,53 1,73 64,7 64,7 64,7 64,7 41,0 53,0 46,9 73,6 81,1 61,4 50,24 77,1 64,85 61,00 67,15 67,13 66,35 63,23 66,50 61,35 64,00 62,50 61,90 59,70 Точки симплекса Худшая точка 1,2,3,4,5 т. 3 1,2,4,5,6 1,2,5,6,7 1,2,5,6,8 1,2,6,8,9 2,6,8,9,10 2,6,8,10,11 2,8,10,11,12 т. 4 т. 7 т. 5 т. 1 т. 9 т. 6 Координаты (условия) 6-й точки симплекса: x16 = 20,3 – 0,3 = 0,3; x26 = 264,39 – 42,6 = 86,2; x36 = 21,16 – 1,36 = 0,96; х46 = 258,8 – 64,75 = 52,9. Записываем условия шестого опыта в табл. 3.8. Проводим опыт в т.6. В симплексе 1,2,4,5,6 выбираем наихудшую точку. Это точка 4. Ее также зеркально отображаем. Подобная процедура повторяется до тех пор, пока не достигнем оптимального результата (т. 12). Вопросы для самоконтроля 1. 2. 3. В чем заключается суть симплексного метода планирования и оптимизации? В чем преимущество симплекс–планирования? Каким образом можно определить, что пришли в оптимальную область? 111 4. МЕТОДЫ ОПТИМИЗАЦИИ В ХИМИЧЕСКОЙ ТЕХНОЛОГИИ 4.1. Основные понятия и определения Конечной целью моделирования химико-технологического процесса (ХТП) является его лучшая реализация или его оптимизация [1, 9]. Оптимизация – это целенаправленная деятельность, которая заключается в получении наилучших результатов (значений параметров объектов) при соответствующих условиях. Оптимизация заключается в нахождении экстремума рассматриваемой функции или оптимальных условий проведения технологического процесса. Для оценки оптимума необходимо прежде всего выбрать критерий оптимизации. Критерием оптимизации (оптимальности) называется количественная оценка оптимизируемого качества объекта. Это главный признак эффективности решения оптимизационной задачи. В зависимости от конкретных условий в качестве критерия оптимальности можно выбрать технологический критерий (например, максимальный выход продукции с единицы объема аппарата), а также экономический критерий (например, минимальная стоимость продукта при заданной производительности). Требования к критерию оптимальности 1. Критерий оптимальности должен быть единственным. 2. Критерий оптимальности должен выражаться числом. На основании выбранного критерия оптимальности составляется целевая функция (функция выгоды), которая представляет собой зависимость критерия оптимальности от параметров, влияющих на его значение. Целевая функция – это критерий оптимальности, рассматриваемый как функция входных параметров: F = F(x1,x2,….,xn). (4.1) Чем больше или меньше F, тем лучше. Следовательно, оптимум – это экстремум (max или min) целевой функции, а задача оптимизации сводится к нахождению экстремума. Оптимизирующие параметры – это те входные параметры системы, которые в процессе оптимизации относят к управляющим и которые применяются для оптимизации процесса. 112 Ограничения – это условия, которые необходимо соблюдать независимо от того, как их соблюдение повлияет на величину критерия оптимальности. Примеры возможных ограничений:  по количеству и качеству сырья и продукции;  по условиям технологии: а) например, в качестве управляющего параметра выбрана температура, которая не может быть выше той, при которой портится (спекается) катализатор; б) не можем менять размер аппарата; в) управляющий параметр – объемная скорость;  величина расхода смеcи ограничивается мощностью насоса;  по экономическим соображениям (капитальные затраты не должны превышать выделенной суммы);  по вопросам охраны труда и окружающей среды. По математическим признакам ограничения разделяют:  на ограничения типа равенств, которые устанавливают определенные значения того или иного фактора: xi = ai (например, задаются значения по составу сырья, размеры аппарата и т. д.);  ограничения типа неравенств, которые определяют пределы изменения параметров процесса. Например, fiai(производительность не ниже заданной); alflbl (температура в определенном интервале); fkbk (температура не выше той, которую выдержит материал). Постановка задачи оптимизации: 1. Необходимо создать математическую модель объекта оптимизации. 2. Выбрать критерий оптимальности, оптимизирующие параметры и сформировать функцию цели. 3. Установить возможные ограничения, которые должны накладываться на переменные. 4. Выбрать метод оптимизации, который позволит найти экстремальное значение искомых величин. Таким образом, математически решить задачу оптимизации – значит определить оптимум функции цели (4.1). Различают задачи статической оптимизации для процессов, протекающих в установившихся режимах, и задачи динамической оптимизации при неустановившихся режимах процесса. 113 4.2. Систематизация методов оптимизации При решении конкретной задачи оптимизации исследователь должен выбрать метод, приводящий к конечным результатам с наименьшим объемом вычислений. Выбор того или иного метода в значительной степени определяется постановкой задачи оптимизации, а также математической моделью объекта оптимизации. Основные методы оптимизации, наиболее широко используемые в химической технологии, можно разделить на несколько групп: 1. Аналитические:  методы исследования функций классического анализа применяются для детерминированных процессов с критерием оптимальности в виде дифференцируемых функций;  метод множителей Лагранжа – для задач с ограничениями типа равенств с критерием оптимальности в виде дифференцируемых функций;  вариационные методы – для задач с критерием оптимальности в виде функционала, расчет оптимальных температурных профилей химических реакторов, оптимальных режимов периодических процессов;  принцип максимума Понтрягина – класс задач с объектами, которые описываются дифференциальными и конечными уравнениями, расчет оптимального управления в задачах регулирования. 2. Методы математического программирования:  геометрическое программирование: процессы с математическим описанием в виде алгебраических функций-полиномов;  динамическое программирование: многостадийные процессы с критерием оптимизации в виде аддитивной функции (секционированные реакторы, каскад аппаратов и т. д.);  линейное программирование: процессы, которые описываются линейными алгебраическими уравнениями с критерием оптимальности в виде линейной функции. 3. Градиентные. Объект оптимизации – сложные процессы химической технологии, отдельные объекты и каскады аппаратов (оптимизация нелинейных и линейных функций с нелинейными и линейными ограничениями). 4. Статистические. Объекты оптимизации не имеют детерменированного описания. При составлении алгоритмов и программ с использованием методов оптимизации целесообразно придерживаться модульного принципа, т. к. требуется многократное обращение к расчету целевой функции. 114 4.3. Статистические методы оптимизации Если в ходе исследования химико-технологического процесса не удается сформировать его детерминированную модель, т. е. неизвестен вид целевой функции, то задача оптимизации решается экспериментально-статистическими методами. Оптимум находят экспериментальным путем. Существует две области изменения выходного параметра y:  область, удаленная от оптимума, в которой происходит значительное изменение выходного параметра y;  почти стационарная область (область, близкая к экстремуму), в которой практически не происходит изменения y. После того как область, удаленная от оптимума, описана линейным уравнением, используем его для оптимизации, т. е. для движения к оптимальной области. Для исследования области, удаленной от оптимума, используются градиентные методы. Рассмотрим метод Бокса-Уилсона. Метод крутого восхождения по поверхности отклика (Бокса-Уилсона) В 1951 г. Д. Бокс и К. Уилсон предложили использовать для поиска оптимальных условий сложных процессов результаты полного или дробного факторного эксперимента [1, 25]. Сформулируем задачу оптимизации. Определить координаты оптимальной (экстремальной) точки опт (x1 , x2опт, …, xnопт) поверхности отклика y = f(x1, …, xn). Метод градиента предусматривает движение к оптимуму по наикратчайшему пути. Метод Бокса-Уилсона – это «шаговый» метод движения по поверхности отклика. Рассмотрим поиск оптимума, когда на процесс влияют два фактора: n = 2 (рис. 4.1). В окрестности точки М (область, удаленная от оптимума) ставится эксперимент по схеме ПФ или ДФ планирования для локального описания поверхности отклика в окрестности т. М, линейным уравнением регрессии yˆ  b  b x  b x2 . 0 1 2 (4.2) Если линейное уравнение адекватно, от центра плана начинают движение к оптимуму по поверхности отклика в направлении градиента: 115 f yˆ   b1; x1 x1 f yˆ   b2 . x2 x2 Движение к экстремуму продолжают до тех пор, пока наблюдается увеличение (поиск максимума) параметра оптимизации y. Движение к экстремуму прекращают в следующих случаях:  если значения факторов или функции отклика выходят за пределы допустимых значений;  если достигнут экстремум критерия оптимальности y. В первом случае поиск оптимума прекращают. Во втором – переносят центр планирования в точку С, до которой дошли по градиенту, и в области локального экстремума функции y выполняют эксперимент по схеме ПФЭ или ДФЭ для математического описания поверхности отклика в окрестности т. С. Если удается вновь получить адекватное математическое описание функции (4.2), то продолжают оптимизацию методом крутого восхождения. Если в области оптимума (т. С) не удалось получить адекватного линейного уравнения регрессии (4.2), то переходят к планированию эксперимента второго порядка для получения математического описания функции цели полиномом второй степени. Полученная поверхность исследуется для локализации экстремума. Метод Бокса-Уилсона используется только для одной экстремальной функции. Рассмотрим расчет значений шагов движения к оптимуму для каждого фактора. При постановке опытов величина шагов должна быть пропорциональна произведению коэффициентов bi на интервал варьирования фактора. Расчет «шагов» при движении по градиенту проводят следующим образом: 1. Один из факторов выбирают за базовый. Вычисляют произведение коэффициента регрессии на соответствующий интервал варьирования, например b1x1: а = (b1x1). 2. Для базового фактора выбирают шаг крутого восхождения ha, оставляя старый интервал варьирования xi либо вводя новый – более мелкий. Обычно ha ≤x1. 3. Производят расчет шага для каждого фактора по уравнению b x hi  i i  ha , a 116 где коэффициенты берутся со своими знаками. Таким образом, знак шага по каждому фактору совпадает со знаком соответствующего коэффициента регрессии. Рассчитанные шаги по каждому целесообразно округлить. Движение начинают от основного уровня. При первом шаге, т. е. опыте, факторы получают значения, равные основному уровню, плюс рассчитанные шаги варьирования: xi = x0 + hi. На каждом последующем шаге значения факторов изменяют на величину шага варьирования. Движение продолжают до тех пор, пока наблюдается увеличение (или уменьшение) функции отклика. Пример. Методом крутого восхождения получить максимальный выход продукта химической реакции. В качестве независимых параметров выбраны концентрация – x1 и температура – x2. Условия эксперимента приведены в табл. 4.1. Таблица 4.1 Фактор x1 (% масс) x2 (°С) Уровни Верхний Нижний Интервалы варьирования xi Основные уровни xi0 40 20 10 30 80 40 20 60 Был проведен эксперимент по схеме ПФЭ (два параллельных опыта) для описания поверхности отклика линейным уравнением y  b  b x  b x . 1 1 2 2 Матрица планирования и результаты опытов приведены в табл. 4.2. Таблица 4.2 N x1 x2 yэ 1 2 3 4 + + – – + – + – 38,4 42,6 43,3 32,2 На основании результатов эксперимента рассчитаны коэффициенты регрессии по формулам (3.35), получено линейное уравнение регрессии следующего вида: y  36,6  5,56 x  1,175 x . (4.3) 1 117 2 Был выполнен статистический анализ по методике (3.36) – (3.42), изложенной в разд. 3.3.1.2. Проверка дисперсии на однородность по критерию Кохрена (3.38) и оценка значимости коэффициентов регрессии по критерию Стьюдента (3.40) показали, что дисперсия однородна и все коэффициенты значимы. Проверка уравнения (4.3) на адекватность, выполненная по критерию Фишера (3.42), показала, что уравнение адекватно и, следовательно, может быть использовано для крутого восхождения по поверхности отклика. Выполним расчет шагов крутого восхождения для каждого фактора (табл. 4.3). Таблица 4.3 Фактор x1 (% масс) x2 (°С) Основные уровни xj0 Интервал варьирования xi Коэффициент bi bixi Шаг К bixi 30 10 5,56 55,6 5,56 60 20 –1,175 –23,5 –2,35 * К – коэффициент (в данном случае 0,1). Движение начинаем от основного уровня (30 %, 60 °С), и на первом шаге факторы принимают значения xi1 = xi0 + hi. Результаты экспериментов приведены в табл. 4.4. Таблица 4.4 Результаты экспериментов по методу крутого восхождения N 1 2 3 4 5 6 1 x1 35,56 41,12 46,68 52,24 57,80 63,36 35,56 x2 57,65 55,30 52,95 50,60 48,25 45,90 57,65 y 40,40 43,56 48,73 52,72 58,00 50,62 40,40 Видим, что в пятом опыте получен максимальный выход продукта реакции. В шестом опыте наблюдается убывание параметра оптимизации y. Переносим центр плана в лучшую точку (т. 5 – т. С на рис. 3.6) и реализуем ПФЭ для двух факторов. Определяем новые условия опытов (табл. 4.5): 118 Таблица 4.5 Фактор x1 (% масс) x2 (°С) Уровни верхний нижний Основной уровни Интервал варьирования xi 57,8 10 67,8 47,8 48,25 20 68,25 28,25 Выполняем эксперимент по схеме ПФЭ 22. Матрица планирования и результаты эксперимента приведены в табл. 4.6. Таблица 4.6 N x0 x1 x2 y 1 2 3 4 +1 +1 +1 +1 +1 –1 +1 –1 +1 +1 –1 –1 56,4 56,6 62,7 70,5 В результате регрессионного анализа уравнений (3.36) – (3.42) получаем, что линейное уравнение регрессии неадекватно эксперименту. Следовательно, достигнута область, близкая к экстремуму (почти стационарная область). Если бы уравнение вновь было бы адекватным, то от основного уровня вновь продолжили бы движение по градиенту (но шаги в этом случае делаются меньшими, т. к. при приближении к оптимуму кривизна поверхности возрастает). 4.4. Аналитические методы Аналитические методы являются классическими методами отыскания экстремального значения функции (min или max). Они применяются в основном в тех случаях, когда известен аналитический вид оптимизируемой функции F от независимых переменных xi и число переменных xi невелико. При большом числе переменных возникает так называемый барьер многомерности и применение аналитических методов становится затруднительным [1, 9]. Аналитический поиск экстремального значения целевой функции F(x1, x2, …, xn) сводится к приравниванию нулю её частных производных: dF/dxi = 0 i = 1, 2,…, n. 119 Необходимые и достаточные условия существования экстремума функции одной переменной Необходимые условия существования экстремума непрерывной функции F(x) (при отсутствии ограничений) могут быть получены на основании анализа первой производной dF/dx. Функция F(x) может иметь экстремальное значение в тех точках оси х, где производная dF/dxi равна нулю или не существует. Равенство нулю производной графически означает, что касательная к кривой F(x) в этой точке параллельна оси х (рис. 4.1). Рис. 4.1 Покажем случаи, когда производные в точках экстремума не существуют. F (x) F (x) а б Рис. 4.2. Типы экстремумов 120 В точках min и max существует конечный разрыв производной dF/dx (рис. 4.2, а). На рис. 4.2, б приведен случай, когда в точках экстремума значение производной обращается в бесконечность. Происходит бесконечный разрыв производной, при котором её значение изменяется от + до – в т. x1 и от – до + в т. x2. dF  0 , отсутствие производной) – Перечисленные выше условия ( dx необходимые условия существования экстремума. Но их выполнение не означает наличия экстремума функции в данной точке (рис. 4.3). F (x) F (x) Рис. 4.3. Случаи отсутствия экстремума при выполнении необходимого условия его существования Для того чтобы определить, действительно ли в точке существует экстремум, необходимо провести следующие дополнительные исследования: 1. Сравнение значений функции. Рассчитывают значение функции в точке, подозреваемой на экстремум, и в двух близких точках – слева и справа от неё (xi +  и xi – , где  – малая положительная величина). Если оба значения F(xi + ) и F(xi – ) меньше или больше F(xi), то в точке xi – максимум или минимум соответственно. Если же F(xi) имеет промежуточное значение, то в точке F(xi) нет экстремума. 2. Сравнение знаков производных. dF в точках (xi + ) и (xi – ). Определяются знаки производных dx Если знаки различны, то в точке xi – экстремум [если знак меняется с (+) на (–), то в точке xi – максимум, если с (–) на (+) – минимум]. Если знаки совпадают в точках (xi + ) и (xi – ), то точка xi не является экстремальной. 121 3. Исследование знаков высших производных. Этот способ можно применять, если в точках xi (подозреваемs[ на экстремум) существуют производные высших порядков. Вычисляется d 2F . производная dx 2 d 2F d 2F  0 , то в точке xi – максимум, если  0 , то – минимум. Если dx 2 dx 2 Чаще всего пользуются двумя первыми способами, т. к. последний достаточно громоздок. Экстремумы функций многих переменных Решение задачи оптимизации усложняется, если критерий оптимальности является функцией нескольких независимых переменных. Для непрерывной функции F = F(x1, x2,…, xn), имеющей непрерывные производные первого и второго порядков по всем переменным xi (i = 1,…, n), необходимым условием экстремума в точке xi служит равенство нулю частных производных по всем переменным, т.е. точки, в которых может быть экстремум функции, определяются решением системы уравнений F ( x1 , x2 ,..., xn ) i  1,..., n.  0, (4.4) xi Левые части уравнений есть функции факторов x1,…,xn. Поэтому решение системы (4.4) дает оптимальное значение факторов. Если оптимизируется технологический процесс, то этому решению соответствует оптимальный режим. Рассмотрим частные задачи оптимизации ХТП с использованием математических моделей. 4.4.1. Оптимизация реактора идеального смешения В реакторе идеального смешения протекает реакция Определить оптимальное время пребывания реагентов в реакторе, при котором достигается максимальный выход целевого продукта B (рис. 4.4). 122 C C A0 CC CBmax CB CA Рис. 4.4 Составим математическую модель: dC A 1  C A0  C A  k1C A ; dt  dCB 1  CB0  CB  k1C A  k2CB ; dt  dCC 1  CC0  CC  k2CB . dt  Начальные условия: при t = 0 СА(0) = СА0; СВ(0) = СВ0. В стационарном режиме работы реактора dCB  0; dt при C B0  0 уравнение (4.6) примет следующий вид:       (4.5) (4.6) (4.7) 1  CB  k1C A  k2CB  0;  k1C A    CB  k2CB   ; k C  CB  1 A . 1  k2 (4.8) Приравниваем к нулю уравнение (4.5): C A0  C A  k1C A  0; C A0  C A 1  k1  ; CA  C A0 1  k1 . 123 (4.9) Полученное выражение подставим в (4.8): k1C A0 . CB  1  k2 1  k1  (4.10) Для определения оптимального времени контакта (опт), при котором достигается максимальное значение концентрации СВ, необходимо уравнение (4.10) продифференцировать по  и приравнять производную к нулю: dF k1C A0 (1  k1 ) 1  k2   k1C A0  k1 1  k2   k2 1  k1     0. (4.11) 2 d 1  k1 1  k2   Отсюда выразим время контакта: 1  опт  ; k1  k2 k1C A0  F  CBmax   k1 1  k1  k2  1 k1  k2 (4.12)  k2 1   k1  k2     . 4.4.2. Задача поиска оптимальной температуры обратимой химической реакции Протекает химическая реакция Если химическая реакция протекает без побочных стадий, то в качестве критерия оптимальности может быть выбрана скорость реакции. Целевая функции имеет вид F  W  k01  e  E1 RT  C A  k02  e  E2 RT  CB . (4.13) Установим ограничения и выберем оптимизирующие факторы. Критерий оптимальности F зависит от трех параметров: T, CA и CB. Но СА и СВ не могут быть выбраны в качестве оптимизирующих параметров, т. к. они не являются входами системы, а являются результатами реакции, т. е. для увеличения скорости необходимо иметь как можно больше вещества СА и меньше СВ. Цель же процесса противоположная – увеличить концентрацию вещества В и уменьшить концентрацию вещества А. Поэтому концентрации СА и СВ нельзя считать независимыми факторами. 124 Таким образом, есть лишь один независимый параметр, влияющий на функцию цели F – температура. Поэтому настоящая задача является задачей об оптимальной температуре химической реакции. Однако при различных значениях CA и CB влияние температуры может быть различным. Поэтому ставим задачу следующим образом: найти оптимальную температуру химической реакции при фиксированных значениях CA и CB. Таким образом, концентрации CA и CB выступают как ограничения в виде равенств C A/t 0  C A0 ;   CB /t 0  CB0 . Второе ограничение типа неравенств (обязательное): температура не может превысить некоторого максимального значения Tmax: TTmax. Если реакция необратима, т. е. k2 = 0, то в уравнении остается первый член, который с ростом температуры растет неограниченно. В этом случае оптимум определяется ограничением: Tопт = Tmax; W  k1,0  C A  e  E1 RT  k2,0  CB  e  E2 RT ; dW  0; dt k1,0  C A  e k1,0  C A  e  E1 RT  E1 RT  E2 E1 1 E 1   2  k2,0  CB  e RT  2  2  0; R T R T  E1  k2,0  CB  e  E2 RT E1  E2 E1k1,0  C A e   E1  e RT ; E2 k2,0  CB e RT  E  k  C  E  E2 ln  1 1,0 A   1 ;  E k C  RT  2 2,0 B  E1  E2 T .  E1  k1,0  C A  k  ln    E2  k2,0  CB  125  E2 RT  E2 ; 4.5. Численные методы решения оптимизационных задач без ограничений При описании численных методов оптимизации можно выделить два случая поиска оптимума:  одномерный поиск, когда функция цели F зависит от одного оптимизирующего фактора F = F(x);  многомерный поиск, когда оптимизирующих факторов больше одного, F = F(x1,x2,…,xn). 4.5.1. Одномерная оптимизация Оптимизация функции одной переменной – это наиболее простой тип оптимизационных задач. Тем не менее задачи одномерной оптимизации достаточно часто встречаются в практической деятельности инженера [9, 27]. Следует отметить, что все методы одномерного поиска базируются на последовательном уменьшении интервала, содержащего точку экстремума. Об эффективности алгоритмов различных методов можно судить по числу вычислений функции, необходимому для достижения заданной точности. Существует достаточно большое количество методов однопараметрической оптимизации. Рассмотрим некоторые наиболее часто применяемые методы. 4.5.1.1. Метод дихотомии Простейшим методом нахождения экстремума функции одной переменной является метод дихотомии (метод деления отрезка пополам). Пусть функция F (x) унимодальна на отрезке [a, b]. Необходимо найти оптимум функции на этом отрезке (рис. 4.5) с заданной степенью погрешности . Рассмотрим последовательность поиска максимума на отрезке [a, b]. Делим отрезок [a, b] пополам (т. l). Произвольно выбираем малое приращение x () и откладываем его слева или справа относительно т. l: x1 = l + ; x2 = l – , где величина  ≤  (например,  = /2). 126  Рис. 4.5. Поиск оптимума методом половинного деления Рассчитываем значения функции в двух новых точках – F(x1) и F(x2) и сравниваем их. Если F(x1) F(x2), то выбрали бы отрезок [а, x2]. После выбора той или иной половины отрезка задача возвращается к исходным позициям. Опять задан отрезок а, b, на котором надо найти максимум. Поэтому проводим следующий цикл расчета, подобный предыдущему. Процедура вычислений повторяется, пока не выполнится условие b – a  . 4.5.1.2. Метод золотого сечения Одним из наиболее эффективных методов оптимизации, в которых при ограниченном количестве вычислений F(x) достигается наилучшая точность, является метод золотого сечения. В основе данного метода лежит правило геометрического соотношения: отношение длины всего отрезка к большей его части равно отношению большей его части к меньшей. Разделим отрезок l на две части m и l – m, где m, l – m – меньшая и большая части отрезка соответственно: 127  m  ; m m m  (   m ) 2 ; m   2  2m  m 2  0; m2  3m   2  0. Решаем квадратное уравнение относительно m: 2 3 3 9 2  4 2  3  2 ; m       2 2 2  2   5 2 3   5  3  5 ;   2 2 2 m   3 5 ; 2 m = 0,382l = (1 – 0,618)l; l – m = 0,618l. (4.14) Золотое сечение Рассмотрим поиск максимума на отрезке [a, b] (рис. 4.6). Начинаем с деления отрезка [a, b] слева и справа в отношении золотого сечения, получаем точки x1 и x2: x1 = a + 0,382(b – a); x2 = b – 0,382(b – a). В этих точках вычисляются значения функции F(x1) и F(x2) и определяется новый интервал, на котором локализован экстремум. Рис. 4.6. Поиск оптимума методом золотого сечения 128 Согласно рис. 4.6 на отрезке [a, x1] максимума быть не может (если функция унимодальна), поэтому эту часть отрезка отбрасываем, переносим т. a в т. x1 и рассматриваем новый отрезок [a, b]. На этом отрезке уже есть точка (x2), в которой рассчитано значение F(x2). Точка x2 отсекла от прежнего отрезка справа 38,2 %, отсекает от нового (меньшего) 61,8 %. Таким образом, и на новом отрезке т. x2 является точкой золотого сечения. Теперь ее можно назвать точкой x1 и добавить на уменьшенном отрезке только одну – т. x2. Данная процедура продолжается до достижения заданной степени точности: b  a  . Таким образом, на каждом этапе расчета, кроме первого, необходимо рассчитывать значение функции F только в одной точке, что повышает эффективность метода. 4.5.1.2. Метод сканирования Метод сканирования – метод определения оптимального значения функции F, который заключается в сужении интервала [a, b] до заданных размеров. Например, ставится задача поиска оптимальной температуры химической реакции в интервале 100…200 °С (100 °С ≤ Топт ≤ 200 °С). В данном случае необходимо решать задачу оптимизации. Если же интервал значительно меньше (145 °С ≤ Топт ≤ 146 °С), то можно считать, что задача оптимизации решена. Применяется данный метод к непрерывным функциям. Сканированием можно исследовать как функцию одной, так и нескольких переменных. Рассмотрим одномерное сканирование. Пусть на отрезке [a, b] (интервал неопределенности) требуется отыскать экстремум (максимум) целевой функции (рис. 4.7). F(x) Fmax хопт а в 2 x Рис. 4.7. Поиск оптимума методом сканирования 129 Задача поиска экстремума сводится к сужению интервала неопределенности. Зададим количество узлов k на интервале [a, b], в которых будет рассчитываться целевая функция. Рассчитаем величину шага сканирования по выражению ∆x = (b – a)/(k – 1). Интервал поиска [a, b] разбивается на равные участки с шагом ∆x, и во всех точках разбиения определяются значения функции F(xi). Выбирается наибольшее из полученных значений функции (в случае поиска максимума, рис. 4.8). Далее происходит сужение интервала: а = хопт – ∆x; b = хопт + ∆x. Проверяется условие b  a  . Если условие выполняется, то сканирование прекращается; если же не выполняется, то процедура сканирования продолжается на новом, суженном интервале, величина которого составит 2x (новый отрезок [a, b]). 4.5.2. Многомерный поиск оптимума При оптимизации технологических процессов необходимость многомерного поиска оптимума возникает достаточно часто. На основе входных параметров формируется критерий оптимальности и выбирается метод многомерного поиска оптимума целевой функции: F = F(x1,x2,…,xn). (4.15) Рассмотрим один из наиболее часто применяемых методов многомерной оптимизации: метод покоординатного спуска. Покоординатный спуск Рассмотрим поиск минимума (рис. 4.8) целевой функции для случая с двумя факторами: F = F(x1, x2). (4.16) В качестве начального приближения в двумерном пространстве выбираются координаты начальной точки поиска: x10 и x2 0 , т. е. те значения, от которых начинается поиск оптимума. Выбираются величины шагов движения по параметрам: h1 и h2 и малые значения: 1 и 2. Выбор этих величин определяется физическим смыслом задачи. 130 Рис. 4.8. Схема движения к оптимуму методом покоординатного спуска:  – шаг в нужном направлении;  – неудачный шаг Рассчитываем значение функции в начальной точке F(x10, x20) – точка 1. Зафиксируем координату х2 = const и начинаем движение вдоль оси х1 с шагом h1 в сторону уменьшения функции цели (4.16). Движение продолжается до тех пор, пока наблюдается уменьшение функции (на рис. 4.8 – это точки 2, 3, …, 7): Fi + 1 < Fi. В точке 8 получаем значение F большее, чем в точке 7. Поэтому возвращаемся в точку 6, фиксируем координату x1 (х17 = const) и движемся вдоль оси x2 с шагом h2 до тех пор, пока уменьшается функция цели: Fi + 1 < Fi. После очередного неудачного шага (Fi + 1 > Fi) меняем координату, по которой идет поиск минимума функции, и вновь продолжаем движение. Процесс последовательно продолжается до тех пор, пока не будет достигнута заданная точность локализации экстремума, т. е. если шаг по каждому параметру приводит к возрастанию функции цели (поиск минимума), а величина шага меньше или равна заданной степени точности, то расчет прекращается: h1 1; h22. В том случае, если эти условия не выполняются, движение продолжается из лучшей точки с уменьшенной величиной шагов. 131 Вопросы для самоконтроля 1. 2. 3. 4. 5. 6. 7. 8. 9. Дайте понятие оптимизации, критерия оптимальности, оптимизирующих параметров, ограничений. Что включает постановка задачи оптимизации? В каких случаях применяют статистические методы оптимизации? Стратегия метода крутого восхождения по поверхности отклика. В каких случаях для поиска оптимума применяются аналитические методы? Привести этапы оптимизации РИС? Перечислить численные методы оптимизации. Назовите методы одномерной оптимизации. Поиск оптимума методом покоординатного спуска. 132 СПИСОК ЛИТЕРАТУРЫ 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. Кафаров В.В. Методы кибернетики в химии и химической технологии. – М.: Химия, 1985. – 489 с. Холоднов В.А., Дьяконов В.П. Математическое моделирование и оптимизация химико-технологических процессов. Практическое руководство – СПб.: АНО НПО «Профессионал», 2003. – 480 с. Закгейм А.Ю. Введение в моделирование химико-технологических процессов. – М.: Химия, 1982 с. Бондарь А.Г. Математическое моделирование в химической технологии. – М.: Высш. шк., 1973. – 280 с. Гумеров А.Н., Валеев А.Н и др. Математическое моделирование химико-технологических процессов: учеб. пособие. – М.: КолосС, 2008. – 160 с. Пахомов А.Н. Коновалов А.Н. и др. Основы моделирования химико-технологических систем: учеб. пособие. – Тамбов: Изд-во Тамб. гос. техн. ун-та, 2008. – 80 с. Касаткин А.Н. Основные процессы и аппараты химической технологии: учебник для вузов. – М.: АльянС, 2004. – 750 с. Кафаров В.В., Глебов М.Б. Математическое моделирование основных процессов химических производств. – М.: Высш. шк., 1991. – 400 с. Гартман Т.Н., Клушин Д.В. Основы компьютерного моделирования химико-технологических процессов: учеб. пособие для вузов. – М.: ИКЦ «Академкнига», 2006. – 416 с. Расчеты основных процессов и аппаратов нефтепереработки: справочник / под ред. Е.Н. Судакова. – М.: Химия, 1979. – 568 с. Кравцов А.В., Мойзес О.Е., Кузьменко Е.А. и др. Информатика и вычислительная математика: учеб. пособие для вузов. – Томск: ТПУ, 2003. – 243 с. Бенедек, П., Ласло А Научные основы химической технологии: перевод с нем. / под ред. П.Г. Романкова, М.И. Курочкиной. – Л.: Химия, 1970. – 376 с. Плановский А.Н. Процессы и аппараты химической и нефтехимической технологии: учебник / А.Н. Плановский, П.И. Николаев. – 3-е изд., испр. и доп. – М.: Химия, 1987. – 496 с. Хала Э. Равновесие между жидкостью и паром. – M.: Изд. иностр. лит., 1962 133 15. Рид Р., Праусниц Дж., Шервуд Т. Свойства жидкостей и газов. – Л.: Химия, 1982. – С. 110–111. 16. Вольф А.В., Самборская М.А. Технологическое проектирование тарельчатых колонн фракционирования нефти: учеб. пособие – Томск: Изд-во ТПУ, 2009. – 16 с. 17. Танатаров М.А. и др. Технологические расчеты установок переработки нефти. – М.: Химия, 1987. – 350 с. 18. Кафаров В.В. Разделение многокомпонентных систем в химической технологии. Методы расчета. – М.: Московский химикотехнологический институт, 1987. – 84 с. 19. Жоров Ю.М. Моделирование физико-химических процессов нефтепереработки и нефтехимии / Ю.М. Жоров. – М.: Химия, 1978. – 376 с. 20. Яблонский Г.С., Быков В.И., Горбань А.И. Кинетические модели каталитических реакций. – Новосибирск: Наука, 1983. – 254 с. 21. Бесков В.С., Флор К.В. Моделирование каталитических процессов и реакторов. – М.: Химия, 1991. – 252 с. 22. Панченков Г.М., Лебедев В.П. Химическая кинетика и катализ. – М.: Химия, 1985. – 589 с. 23. Киперман С.Л. Основы химической кинетики в гетерогенном катализе. – М.: Химия, 1979. – 349 с. 24. Безденежных А.А. Инженерные методы составления уравнений скоростей реакций и расчета кинетических констант. – Л.: Химия, 1973. – 256 с. 25. Ахназарова С.Л., Кафаров В.В. Оптимизация эксперимента в химии и химической технологии. – М.: Высш. шк., 1978. – 319 с. 26. Налимов В.В., Чернова Н.А. Статистические методы планирования экстремального эксперимента. – М.: Наука, 1965. – 340 с. 27. Кравцов А.В., Новиков А.А., Коваль П.И. Компьютерный анализ технологических процессов. – Новосибирск: Наука, 1998. – 212 с. 134 Учебное издание УШЕВА Наталья Викторовна МОЙЗЕС Ольга Ефимовна МИТЯНИНА Ольга Евгеньевна КУЗЬМЕНКО Елена Анатольевна МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ХИМИКО-ТЕХНОЛОГИЧЕСКИХ ПРОЦЕССОВ Учебное пособие Научный редактор доктор технических наук, профессор А.В. Аристов Корректура Н.Т. Синельникова Компьютерная верстка К.С. Чечельницкая Дизайн обложки Т.А. Фатеева Подписано к печати 20.05.2014. Формат 60×84/16. Бумага «Снегурочка». Печать XEROX. Усл. печ. л. 7,85. Уч.-изд. л. 7,1. Заказ 396-14. Тираж 100 экз. Национальный исследовательский Томский политехнический университет Система менеджмента качества Издательства Томского политехнического университета сертифицирована в соответствии с требованиями ISO 9001:2008 . 634050, г. Томск, пр. Ленина, 30 Тел./факс: 8(3822)56-35-35, www.tpu.ru