close

Вход

Забыли?

вход по аккаунту

?

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

код для вставкиСкачать
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Федеральное агентство по образованию
Государственное образовательное учреждение
высшего профессионального образования
«Казанский государственный технологический
университет»
А.В. Клинов, А.Г. Мухаметзянова
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
ХИМИКО-ТЕХНОЛОГИЧЕСКИХ ПРОЦЕССОВ
Учебное пособие
Казань
2009
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
УДК 66.011+51.74+519.2
ББК 35:22.1
Математическое
моделирование
химикотехнологических процессов: учебное пособие / А.В. Клинов,
А.Г. Мухаметзянова. − Казань: Изд-во Казан. гос. технол. ун-та,
2009. − 144 с.
ISBN 978-5-7882-0774-2
Изложена теория математического моделирования
химико-технологических процессов, рассмотрены типовые
математические модели, применяющиеся при описании химикотехнологических процессов, этапы их построения. Подробно
разбираются методы и модели для определения физикохимических свойств газовых и жидких смесей.
Предназначено для студентов четвертых курсов,
обучающихся по специальности 240802 – «Основные процессы
химических производств и химическая кибернетика», а также
осваивающих другие инженерные специальности. Может
использоваться аспирантами и преподавателями технических
вузов, научными работниками.
Подготовлено на кафедре «Процессы и аппараты
химической технологии».
Печатается по
совета
Казанского
университета.
решению редакционно-издательского
государственного
технологического
Рецензенты: д-р техн. наук, проф. А.Г. Лаптев
д-р техн. наук, проф. Р.И. Ибятов
© Клинов А.В., Мухаметзянова А.Г., 2009
©
Казанский
технологический
университет, 2009
государственный
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ОГЛАВЛЕНИЕ
Введение
1. Основные понятия и принципы моделирования
2. Эмпирический метод построения математических
моделей
2.1. Формулирование цели, выбор факторов и
переменных состояния объекта исследования, виды
уравнений регрессии
2.2. Планирование и проведение экспериментов
2.3. Стaтистическaя обрaботкa результатов и поиск
наилучшей формы aппроксимaции полученных
данных
3. Теоретический метод построения математических
моделей
3.1. Теоретические основы
3.1.1. Законы сохранения
3.1.2. Законы термодинамики и
термодинамического равновесия
3.1.3. Законы переноса и химической кинетики
3.2. Исчерпывающее описание процессов химической
технологии и типовые модели структуры потока
3.2.1. Импульсный ввод индикатора для
определения параметров типовых и комбинированных
моделей структуры потоков
3.2.2. Моделирование теплообменных процессов
3.2.3. Моделирование массообменных процессов
3.2.4. Математические модели химических
реакторов
3.2.5. Влияние структуры потоков на протекание
химических реакций
3.2.6. Влияние переноса массы на протекание
химической реакции и селективность
4. Методы и модели определения физико-химических
свойств газовых и жидких смесей
4.1. Уравнения состояния
4.2. Расчет термодинамических свойств на основе
избыточных функций
4.3. Молекулярно−статистические методы описания
3
5
6
10
12
16
18
29
30
30
34
51
55
62
76
84
99
102
109
115
115
121
126
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
макросвойств газов и жидкостей
Заключение
Библиографический список
4
142
143
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ВВЕДЕНИЕ
В настоящее время в связи с массовым использованием
компьютерной
техники
математическое
моделирование
приобрело повсеместное признание и распространение.
В учебном пособии рассмотрены основные подходы к
математическому моделированию химико-технологических
процессов. В первой главе даны основные понятия метода
математического
моделирования
химико-технологических
процессов, приведена классификация математических моделей.
Вторая глава подробно описывает эмпирический метод
построения математических моделей, который базируется на
статистической обработке результатов экспериментов и поиске
наилучшей формы аппроксимации полученных данных.
Теоретический метод построения математических моделей,
приведенный в третьей главе, подразумевает знание
фундаментальных законов, составляющих теоретическую
основу всех технологических процессов. Поскольку законы
сохранения и переноса подробно рассматривались при изучении
дисциплин «Явления переноса» и «Процессы и аппараты
химической технологии», в настоящем пособии приводятся
только
общие
сведения.
Рассматриваются
законы
термодинамики и их важные следствия, необходимые для
анализа протекания процессов химической технологии, типовые
модели структуры потоков и способы их получения. Приводятся
ряд математических моделей процессов переноса тепла, массы и
химических превращений в аппаратах с различной структурой
потока и модели, учитывающие взаимное влияние этих
процессов друг на друга. Сделан акцент на математические
модели с распределенными параметрами, которые описываются
дифференциальными уравнениями, и граничные условия к ним.
Четвертая глава посвящена методам и моделям
определения физико-химических свойств газовых и жидких
смесей.
В отличие от курса ПАХТ, в учебном пособии не
рассматриваются теории пограничных слоев и межфазного
переноса,
не
обсуждаются
некоторые
особенности
многокомпонентного массопереноса.
5
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
1. ОСНОВНЫЕ ПОНЯТИЯ И ПРИНЦИПЫ
МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ
В настоящее время нельзя назвать область человеческой
деятельности, в которой в той или иной степени не
использовались бы методы моделирования. В первую очередь
это относится к научным исследованиям и производственным
процессам.
Метод моделирования широко применяют при
проектировании,
создании,
внедрении,
эксплуатации
технологических систем, а также на различных уровнях их
изучения, начиная от анализа элементарных явлений и
заканчивая исследованием системы в целом при взаимодействии
с окружающей средой.
Введем основные понятия метода математического
моделирования.
Модель — условный образ объекта исследования,
конструируемый исследователем так, чтобы отобразить
основные характеристики и существенные особенности его
поведения. Такое описание особенно полезно в случаях, когда
исследование самого объекта затруднено или физически
невозможно. Зачастую в качестве модели выступает другой
материальный или мысленно представляемый объект,
замещающий в процессе исследования объект-оригинал. Таким
образом, модель выступает своеобразным инструментом для
познания, который исследователь ставит между собой и
объектом и с помощью которого изучает интересующий его
объект.
В силу многозначности понятия «модель» в науке и
технике не существует единой классификации видов моделей:
классификацию можно проводить по характеру моделей, по
характеру моделируемых объектов, по сферам приложения
моделирования (в технике, физических науках, кибернетике и
т.д.). Тем не менее, применительно к естественным и
техническим наукам принято различать следующие виды
моделей:
• концептуальные модели − совокупность уже
известных
фактов
или
представлений
относительно
исследуемого объекта или системы истолковывается с помощью
6
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
некоторых специальных знаков, символов, операций над ними
или с помощью естественного или искусственного языков;
• физические модели − модель и моделируемый объект
представляют собой реальные объекты или процессы единой
или различной физической природы, причем между процессами
в объекте-оригинале и в модели выполняются некоторые
соотношения подобия, вытекающие из схожести физических
явлений;
• структурно-функциональные модели − моделями
являются схемы (блок-схемы), графики, чертежи, диаграммы,
таблицы, рисунки, дополненные специальными правилами их
объединения и преобразования;
• математические (логико-математические) модели −
построение модели осуществляется средствами математики и
логики;
• имитационные (программные) модели − логикоматематическая модель исследуемого объекта представляет
собой алгоритм функционирования объекта, реализованный в
виде программного комплекса для компьютера.
Разумеется, перечисленные выше виды моделей не
являются взаимоисключающими и могут применяться при
исследовании сложных объектов либо одновременно, либо в
некоторой комбинации.
Моделирование — это процесс построения моделей
(математических или физических) и изучение на их основе
реально существующих предметов, процессов или явлений с
целью получения объяснений этих явлений, а также для
предсказания явлений, интересующих исследователя.
Перенос результатов, полученных в ходе построения и
исследования модели, на оригинал основан на том, что модель в
определенном смысле отображает (воспроизводит, моделирует,
описывает, имитирует) некоторые интересующие исследователя
черты объекта. Моделирование как форма отражения
действительности широко распространенно.
В настоящем учебном пособии будет рассматриваться
математическое моделирование, которое является эффективным
инструментом исследования, прогнозирования и управления
сложными
нестационарными
химико-технологическими
процессами (ХТП). Современная компьютерная техника,
7
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
быстрое развитие вычислительной математики, повсеместное
использование
вычислительной
техники
чрезвычайно
расширило возможности математического моделирования.
Математическая модель, как правило, учитывает лишь те
свойства, которые отражают и определяют основные
особенности
поведения
объекта-оригинала,
а
также
представляют интерес с точки зрения целей и задач конкретного
исследования. Следовательно, в зависимости от целей
моделирования при рассмотрении одного и того же объектаоригинала с различных точек зрения и в различных аспектах
объект-оригинал может иметь различные математические
описания и, как следствие, быть представлен различными
математическими моделями.
Дадим наиболее общие определения математической
модели, математического моделирования, приведем основные
этапы математического моделирования.
Математическая модель − система математических
выражений,
описывающих
характеристики
объекта
моделирования.
Моделирование в химической технологии – это метод
исследования химико-технологических процессов или систем
при помощи построения и изучения их моделей, которые
отличаются от объектов моделирования масштабами или
физической природой происходящих в них явлений, но
достаточно точно (адекватно) отображающих представляющие
интерес свойства этих объектов. Моделирование используют
для решения различных задач, важнейшие из которых: 1)
исследование новых процессов; 2) проектирование производств;
3) оптимизация отдельных аппаратов и технологических схем;
4) выявление резервов мощности и отыскание наиболее
эффективных
способов
модернизации
действующих
производств; 5) оптимальное планирование производств; 6)
разработка
автоматизированных
систем
управления
проектируемыми
производствами;
7)
построение
автоматизированных систем научных исследований.
Математическое моделирование − метод исследования
процессов или явлений при помощи построения их
математических моделей.
8
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Поскольку объектом моделирования здесь будут
выступать ХТП, дадим определение этим процессам.
Химико-технологические
процессы
(ХТП)
–
технологические процессы, связанные с физико-химической и
химической переработкой реагентов в конечные продукты. ХТП
– это сложная физико-химическая система, которая имеет
двойственную
детерминированную
и
стохастическую
(случайную) структуру, переменную, как во времени, так и в
пространстве.
Основными элементами ХТП являются следующие
«элементарные» процессы:
• гидромеханические;
• тепловые;
• массообменные;
• механические;
• химические.
Все эти процессы взаимосвязаны, то есть влияют друг на
друга, что создает известные трудности математического
моделирования ХТП.
Моделирование как метод исследования ХТП включает в
себя следующие основные этапы:
• постановка цели моделирования;
• построение модели;
• идентификация модели;
• проверка адекватности модели и внесение корректив;
• использование ее для исследования свойств и
поведения объекта.
Преимущества мaтемaтического моделирования:
- позволяет осуществить с помощью одного устройства (ЭВМ)
решение целого клaссa зaдaч, имеющих одинaковое
мaтемaтическое описaние;
- обеспечивaет простоту перехода от одной зaдaчи к другой,
позволяет вводить переменные пaрaметры, возмущения и
рaзличные условия однозначности;
- дает возможность проводить моделирование по частям
("элементарным процессам"), что особенно существенно при
исследовании сложных объектов химической технологии;
- экономичнее метода физического моделирования.
9
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Клaссификaция математических моделей
По отношению ко времени различают статические и
динамические модели. Первые инвариантны ко времени, а
вторые являются функцией времени.
По отношению к пространству различают модели с
сосредоточенными
параметрами
и
распределенными
параметрами. Для моделей с сосредоточенными параметрами
характерно постоянство переменных в пространстве. Тогда
математическое описание включает алгебраические уравнения
либо дифференциальные уравнения первого порядка для
нестационарных процессов. Если основные переменные
процесса изменяются как во времени, так и в пространстве или
если указанные изменения происходят только в пространстве, то
модели, описывающие такие процессы, называются моделями с
распределенными параметрами. Их математическое описание
включает обычно дифференциальные уравнения с частными
производными,
либо обыкновенные
дифференциальные
уравнения.
Основные методы построения математических
моделей
Существует два метода составления мaтемaтического
описания:
1) эмпирический (экспериментaльно-стaтистический,
метод «черного ящика»);
2) теоретический (структурный).
На рис.1.1. представлена общая блок-схема построения
математической
модели
некоторого
объекта.
Если
математическая модель строится эмпирическим методом, тогда
в этой блок-схеме должен отсутствовать блок построения
мысленной модели. То есть объект представляется в виде
«черного ящика», природа и сущность его не известны, модель,
строящаяся
на
основании
статистической
обработки
экспериментальных данных, будет описана в следующих
разделах.
10
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Рис. 1.1. Общая блок-схема построения математической
модели процесса
2. ЭМПИРИЧЕСКИЙ МЕТОД ПОСТРОЕНИЯ
МАТЕМАТИЧЕСКОГО ОПИСАНИЯ
Эмпирический метод в основном используется, когда
процесс мало изучен или ничего неизвестно о его природе. Этот
метод также позволяет получить мaтемaтическое описaние
действующего объекта без исследования и понимания его
внутренней структуры.
При
использовании
эмпирических
методов
мaтемaтическое описaние составляется следующим образом:
1. Формулирование цели, выбор факторов и переменных
состояния объекта исследования, планирование экспериментов.
2. Проведение экспериментов методом «черного ящика»,
то есть изучение реакции объекта на рaзличные возмущения.
11
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
3. Стaтистическaя обрaботкa результатов и поиск
наилучшей формы aппроксимaции полученных данных.
Здесь можно выделить следующие основные этапы:
• проверка опытов на воспроизводимость;
• определение погрешности эксперимента;
• проверка степени корреляции входных и выходных
характеристик;
• построение математической структуры модели
(структурная идентификация);
• определение параметров эмпирической модели
(параметрическая идентификация);
• проверка значимости коэффициентов регрессии;
• проверка на адекватность.
4. Проведение исследований на основе полученной
модели.
2.1. Формулирование цели, выбор факторов и переменных
состояния объекта исследования, виды уравнений регрессии
Внешние связи любой системы можно представить в
виде схемы (рис. 2.1).
Рис. 2.1. Модель «черного ящика»: X,U,Z − вектора
входных параметров;
−
Y вектор выходных параметров
Здесь Х, U, Z – входные факторы:
«Х» – факторы, которые мы можем контролировать и
регулировать (расходы и составы потоков веществ, температура,
давление в потоках, термодинамическое состояние, частота
оборотов мешалки и т.д.);
12
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
«U» – факторы, которые мы можем контролировать, но не
можем регулировать (размеры аппарата, состав сырья);
«Z» – факторы, которые мы не можем ни регулировать, ни
контролировать; вектор неопределенных параметров (часть
входных переменных и внутренних параметров системы,
значения которых мы не знаем точно) и возмущающих
воздействий.
«Y» – выходные факторы, которые называются переменными
состояния, или функцией отклика. Различают экономические и
технологические
переменные
состояния.
В
качестве
экономических используют производительность, себестоимость
и другие показатели. Технологическими переменными являются
качество продукта, выход целевого продукта, надежность
получаемых изделий и др. Объект исследования может иметь
несколько переменных состояния, которые следует сократить до
минимума. Опыт показывает, что в большинстве случаев
удается ограничиться одной переменной состояния, и тогда
вектор Y превращается в скаляр.
Если переменных состояния несколько, то эксперимент
проводится по каждой из них, а затем решается компромиссная
задача.
При выборе переменной состояния необходимо
учитывать следующие требования:
• переменная состояния должна иметь количественную
характеристику, то есть измеряться;
• переменная состояния должна однозначно измерять
эффективность объекта исследования; это требование
эквивалентно корректной постановке задачи;
• переменная состояния должна быть статистически
эффективной, то есть обладать, возможно, меньшей дисперсией
при проведении опытов; это позволяет хорошо различать
опыты.
Правильный выбор переменной состояния объекта
исследования повышает шансы экспериментатора на успех.
К факторам и переменным состояния одновременно
предъявляется ряд требований:
• факторы и переменные состояния должны иметь
области определения, заданные технологическими или
принципиальными ограничениями;
13
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
• между факторами и переменными состояния должно
существовать однозначное соответствие; оно позволит в
основном эксперименте построить математическую модель
объекта исследования и решить поставленную задачу
эксперимента.
Точность измерения и управления входных факторов
должна быть известна и достаточно высока (хотя бы на порядок
выше точности измерения выходной переменной); очевидно,
что низкая точность измерения факторов уменьшает
возможности воспроизведения эксперимента.
Как правило, для каждой отдельной переменной
состояния строится своя модель в зависимости от всех входных
факторов. При составлении математической модели Х и U
считаются заданными или точно определенными величинами. В
отличие от них переменная состояния Y определяется с
некоторой ошибкой и поэтому считается случайной величиной.
В этом случае мaтемaтическое описaние для некоторой
переменной состояния представляет собой уравнение вида
Y = F(U, X, Z)
(2.1)
где Y =
1 n
∑ Yi ; n – число параллельных опытов, в которых
n i =1
значения входных переменных задавались одинаковыми.
Уравнения (2.1) определяют зависимость выхода от всех
входных воздействий. Но установить вид функций F
принципиально невозможно − возмущения Z нам неизвестны.
Однако в большинстве случаев уравнение (2.1) можно
представить в виде
Y = F(U, X) + f (Z).
(2.2)
Здесь функция рaзбитa на два слaгaемых: зависимость F
от контролируемых пaрaметров и погрешность «шум» f. Теперь
зaдaчa ставится таким образом: установить вид функции F и
оценить «шум» f.
Под мaтемaтической моделью будем понимать именно
Y = F(U, X).
(2.3)
14
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Математическая модель процесса всегда содержит
некоторое
количество параметров,
значения
которых
определяются из экспериментальных данных и не зависят от
регулируемых входных факторов. Факторы U, которые
невозможно регулировать или изменять в процессе
исследований, можно также принять за параметр модели или
считать, что их влияние неявным образом будет учитываться
другим параметром модели.
Тогда уравнение (2.3) запишется как
Y = F( X, A),
(2.4)
где A = [ a1 , a 2 ,..., a m ] , m – число параметров эмпирической
модели; X=[x1,x2,…,xk], k − число входных факторов.
Уравнение (2.4) устaнaвливaет связь между выходными и
входными факторами и нaзывaется урaвнением регрессии.
Конкретный вид функциональной зависимости (2.4) и значения
коэффициентов А определяются из опытных данных.
В общем случае различают два вида уравнений
регрессии (эмпирических моделей): нелинейные по параметрам
А, статистический анализ которых осуществляется методом
«нелинейной регрессии», и линейные по параметрам А,
статистический анализ которых проводится методом «линейной
регрессии».
Линейные по параметрам модели могут быть
представлены в следующем общем виде:
m
Y = ∑ a jϕ j (x),
j=1
(2.5)
где ϕ j (x) – линейные или нелинейные функции входных
переменных. В простейшем случае Y = a1 + a 2 x получаем
уравнение линейной регрессии.
Определение параметров (коэффициентов) линейных
моделей и их регрессионный анализ существенно проще, чем
нелинейных моделей, поэтому нелинейные модели, по
возможности, стараются линеаризовать и привести к виду (2.5).
15
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Пользуясь статистическими методами и учитывая
конечность экспериментальных данных, можно получить
оценки коэффициентов регрессии А в уравнении (2.5).
Очевидно, что число экспериментальных данных при этом
должно быть как минимум не меньше числа параметров.
Обычно внaчaле рaссчитывaют более простые
многочлены, отклонение опытных точек от рaсчетных знaчений
срaвнивaют со случaйной ошибкой эксперимента. Если обе
величины
одного
порядка,
то
описaние
считают
удовлетворительным. Если отклонение нельзя объяснить
случaйной ошибкой, то рaссчитывaют более сложный
многочлен. По мере увеличения порядка многочлена точность
описания
возрaстaет,
но
одновременно,
во-первых,
увеличивается требуемое число опытов для нахождения
коэффициентов многочлена, a во-вторых, усложняется
трaктовкa модели.
2.2. Планирование и проведение экспериментов
Урaвнения регрессии можно получить одним из трех
способов:
1. Пассивный эксперимент.
2. Активный эксперимент.
3. Определение реакции объекта на стaндaртные
возмущения.
Пассивный эксперимент − производится сбор и aнaлиз
информации о состоянии технологических пaрaметров объекта
без специального изменения входных пaрaметров процесса. При
пассивном эксперименте исследователь лишь регистрирует
случайные входные воздействия, возникающие при нормальной
эксплуатации объекта, и реакцию объекта на эти воздействия.
Достоинства данного метода:
• Практически полностью отсутствуют зaтрaты на
эксперимент.
Недостатки:
• В нормальных условиях эксплуaтaции колебания
технологического
режима
невелики,
и
поэтому
экспериментальные точки близки друг к другу. В этих условиях
16
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
на точность описания могут сильно повлиять случайные
ошибки.
• Необходимо иметь достаточно большое количество
экспериментальных данных.
Активный эксперимент состоит в целенaпрaвленном
изменении входных пaрaметров технологического процесса. В
основе этого метода лежит плaнировaние эксперимента.
Практически все процессы химической технологии
являются сложными, и на покaзaтели процесса окaзывaют
влияние большое число фaкторов. Возможны два подхода к
исследованию таких многофакторных систем. Первый основан
на том, что исследование объекта рaзбивaется на серии, в
каждой из которых исследуется изменение только одного
пaрaметрa при фиксированных остальных. Второй подход
основан на построении плaнa эксперимента, который
предусмaтривaет изменение всех влияющих фaкторов. Такой
эксперимент нaзывaют многофакторным.
Достоинством первого подхода является его наглядность
и простота интерпретации получаемых результатов. Второй
подход значительно эффективнее − при том же объеме
экспериментальных исследований и той же точности опытов
получается существенно большая точность результатов.
При определении реакции объекта на стaндaртные
возмущения на вход подается кaкой-то стaндaртный сигнал −
единичный импульс, ступенчaтое либо синусоидальное
изменение входного пaрaметрa (рис. 2.2).
Исследование объекта при нанесении стaндaртных
возмущений заметно облегчает обработку получаемой
информации. Этим способом в основном пользуются при
изучении динамики (переходных хaрaктеристик) объекта, при
определении гидродинамической обстановки и др.
17
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Рис. 2.2. Определение реакции объекта на стандартные
возмущения: t − время; Свх − концентрация на входе
ключевого компонента
2.3. Стaтистическaя обрaботкa результатов и поиск
наилучшей формы aппроксимaции полученных данных
Проверка опытов на воспроизводимость
Проводят N серий из n параллельных опытов в
рассматриваемой области изменения влияющих факторов. В
каждом параллельном опыте j-ой серии значения входных
переменных задаются одинаковыми, но для каждой серии они
меняются. Результаты опытов заносятся в табл. 2.1.
Таблица 2.1
Номер
Результаты параллельных опытов
серии
опытов
1
Y11
Y12
..
Y1n
yj,средн.
sj2
s12
2
Y21
Y22
..
Y2n
3
Y31
Y32
..
Y3n
Y1
Y2
Y1
j
Yj1
Yj2
..
Yjn
Yj
sj2
N
YN1
YN2
YNn
YN
s N2
1. Определяется среднее арифметическое значение
функции отклика для любой серии опытов j (j=1,N):
Yj =
(2.6)
18
1 n
∑ Yji .
n i =1
s22
s32
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
2. Оценивается
параллельных опытов:
дисперсия
s 2j =
для
каждой
серии
1 n
(Yji − Y j ) 2 ,
∑
n − 1 i =1
(2.7)
где (n−1) – число степеней свободы − общее число измерений за
вычетом числа оценок, уже рассчитанных по этим измерениям и
применяемых при расчете рассматриваемой характеристики.
3. Определяется расчетное значение критерия Кокрена:
Gp =
max(s 2j )
N
∑s
i =1
.
2
i
(2.8)
4. По специальным таблицам [1] определяют табличное
значение критерия Кокрена − GT. Оно зависит от доверительной
вероятности P (как правило, P = 0,95), от N и от f = n−1.
Если выполняется условие Gp ≤ GT , то опыты считаются
воспроизводимыми.
Определение погрешности эксперимента
Для определения погрешности эксперимента проводится
оценка дисперсии воспроизводимости по формуле
S2Y =
1 N 2
∑ si .
N i =1
(2.9)
Проведение структурно-регрессионного анализа
Методы корреляционного и регрессионного анализов
применяются для выявления и описания зависимостей между
случайными величинами по экспериментальным данным.
Проверяется степень корреляции входных и выходных
характеристик, а также определяется вид функциональной
зависимости Х от Y.
19
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
О наличии или отсутствии корреляции двух случайных
величин качественно можно судить по виду поля корреляции,
нанеся точки (хi,уi) на координатную плоскость. Поля
корреляции случайной величины представлены на рис. 2.3.
В общем виде задача выявления и оценки силы
стохастической связи не решена. Существуют показатели,
оценивающие те или иные стороны стохастической связи.
Основной числовой характеристикой связи двух
случайных величин X и Y является ковариационный
(корреляционный) момент
K xy = (x − x)(y − y).
(2.10)
Рис. 2.3. Поле корреляции случайной величины
20
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
а − сильная отрицательная корреляция; б − сильная
положительная корреляция; в − слабая положительная
корреляция;
г − отсутствие корреляции
•
Назовем случайную величину x = x − x центрированной,
•
ее значение x есть отклонение x от среднего значения x . Тогда
момент K xy равен среднему значению произведения двух
•
•
центрированных случайных величин x, y и характеризует
тесноту взаимной связи X и Y. Действительно, пусть X и Y
независимы, тогда среднее от произведения этих двух величин
будет равно произведению средних:
K xy = (x − x)(y − y) = (x − x)(y − y) = (x − x)(y − y) = 0.
(2.11)
Обратное утверждение – если K xy = 0 , то X и Y
независимы − в общем случае доказать нельзя. Более того,
можно найти пример, в котором K xy = 0 , а X и Y - зависимые
.
случайные величины. Если x мало, то и K xy также будет
малым, хотя X и Y в достаточной мере коррелированны.
Следовательно, величина K xy не является исчерпывающей
характеристикой силы связи X и Y. Оказывается, что K xy
характеризует только линейную зависимость случайных
величин X и Y. Если K xy = 0 , то можно говорить о том, что X и
Y линейно независимы (не коррелированны).
Размерность K xy не имеет физического смысла, что
делает невозможным сравнение ковариационных моментов для
двух пар случайных величин X и Y, поэтому чаще всего
используется нормированная величина K xy - коэффициент
корреляции:
21
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
rxy =
K xy
sxs y
,
(2.12)
где s 2x , s 2y − дисперсии X и Y; s x = s 2x и s y = s 2y .
s 2x =
N
1
(x i − x)2 ,
∑
(N − 1) i =1
s 2y =
N
1
(yi − y) 2 .
∑
(N − 1) i =1
(2.13)
Для
независимых
случайных
величин
rxy = 0 , но
rxy = 0 может быть и для некоторых зависимых величин,
которые при этом называются некоррелированными.
Коэффициент корреляции характеризует не всякую
зависимость, а только линейную. Линейная вероятностная
зависимость случайных величин X и Y заключается в том, что
при возрастании X случайная величина Y имеет тенденцию
возрастать (при rxy > 0 ) или убывать (при rxy < 0 ) по
линейному закону. Величина rxy характеризует степень тесноты
линейной зависимости 1 ≤ rxy ≤ 1 .
При rxy = ±1 величины X и Y связаны между собой
линейной функциональной зависимостью, а при 1 < rxy < 1 связь
между этими величинами стохастическая. При rxy = 0 между X
и Y линейной корреляционной связи не существует, однако
может существовать нелинейная регрессия.
После определения количественной оценки тесноты
связи необходимо решить задачу выборa модели, которая
возникaет при нaличии для одного и того же объекта клaссa
моделей. Задача ставится таким образом: по данной выборке
объема n найти уравнение приближенной регрессии и оценить
допускаемую при этом ошибку. Выбор модели является одним
из вaжнейших этaпов моделировaния. В конечном счете
преимущество той или иной модели определяет критерий
прaктики, понимaемый в широком смысле. При выборе модели
следует исходить из рaзумного компромисa между сложностью
модели, полнотой получaемых с ее помощью хaрaктеристик
22
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
объекта и точностью этих хaрaктеристик. Тaк, если модель
недостaточно точнa, то ее нужно дополнить, уточнить
введением новых фaкторов, может тaкже окaзaться, что
предложеннaя модель слишком сложнa и те же результaты
можно получить с помощью более простой модели.
Расчет коэффициентов регрессии
Коэффициенты
регрессии
рассчитываются
методом
наименьших квадратов. По методу наименьших квадратов
можно обрабатывать любые экспериментальные данные, однако
оптимальность этой процедуры доказывается только для
нормального распределения. Основное условие метода
формулируется следующим образом: коэффициенты регрессии
определяются на основании минимизации суммы квадратов
отклонений между экспериментальными y эj , и рассчитанными
по уравнению регрессии y pj значениями функции отклика:
N
F = ∑ (y эj − y pj )2 → min.
j=1
(2.14)
Если уравнение регрессии содержит m параметров а,
необходимым условием минимума функции F является
выполнение равенств
 ∂F
 ∂a = 0
 1
 ∂F
=0

 ∂a 2
...

 ∂F
 ∂a = 0
 m
.
(2.15)
Для линейных или линеаризованных по параметрам
моделей условие (2.5) приводит к системе линейных уравнений
относительно искомых параметров а. Действительно подставим
23
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
уравнение регрессии вида (2.5) в (2.15), тогда k−ое уравнение
системы (2.15) будет иметь вид
∂F
∂
=
∂a k ∂a k
2
 э m

∑
 y j − ∑ a i ϕi (x j )  = 0.
j=1 
i =1

N
(2.16)
После дифференцирования получим:
N
m


2∑  y эj − ∑ a i ϕi (x j )  ϕk (x j ) = 0.
i =1
j=1 

(2.17)
Перегруппируем члены:
m
N
N
∑ a ∑ ϕ (x )ϕ (x ) = ∑ y ϕ (x ).
i =1
i
i
j=1
j
k
j
j=1
j
k
j
(2.18)
Видно, что это линейное (относительно искомых
параметров) уравнение.
Введем в рассмотрение информационную матрицу I,
элементы которой определяются согласно
N
I ik = ∑ ϕi (x j )ϕk (x j ).
j=1
(2.19)
Эта квадратная симметричная матрица и значения ее
элементов зависят только от входных переменных и
конкретного вида функции ϕ j (x) , причем информационную
матрицу
можно
представить
в
виде
произведения
трансформированной и исходной матрицы входных переменных
Ф:
(2.20)
I = Ф TФ.
Матрица, зависящая от входных переменных, имеет вид
ϕ1 (x1 )...ϕm (x1 ) 
Ф = ..........................  .
ϕ1 (x N )...ϕm (x N ) 
(2.21)
24
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Тогда систему линейных уравнений (2.18) можно
переписать в матричном виде
IA = B ,
(2.22)
где
B = ФТ y э ,
 y1э 
 
y э =  ...  − вектор
 yэ 
 n
значений выходной переменной.
Вектор значений параметров
согласно уравнениям
экспериментальных
модели
определяется
A =I -1B =(Φ T Φ)-1Φy э .
(2.23)
Условия существования единственного решения системы
уравнений (2.23) это
D ( ФТФ )
−1
≠ 0,
(2.24)
Если детерминант будет равен нулю, то система
уравнений имеет бесконечное множество решений. Так же
необходимо помнить, что число искомых параметров должны
быть не больше числа экспериментальных точек (m≥n), иначе
решение получить невозможно.
Если модель нелинейна по параметрам, то такой метод
решения приведет к системе нелинейных уравнений. Решение
подобных систем намного сложнее и может привести к
нескольким наборам параметров.
После
определения
коэффициентов
регрессии
определяют значимость этих коэффициентов. Оценку
значимости коэффициентов модели выполняют по критерию
Стьюдента. Для каждого проверяемого коэффициента
a j вычисляют величину
t j = a j / σa j
25
,
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
где σa j =
 ∂a j  2
SY .
i =1 
i 
m
∑  ∂y
Дисперсия воспроизводимости S2Y определяется по
формуле (2.9).
Если tj больше табличного tT (табличное значение
критерия Стьюдента [1], которое находится по числу степеней
свободы и доверительной вероятности), то коэффициент
aj значительно отличается от нуля.
Незначимые коэффициенты отбрасываются из уравнения
регрессии,
после
чего
оставшиеся
коэффициенты
пересчитываются заново, поскольку они закоррелированы друг
с другом.
Полученное уравнение проверяется на адекватность.
Проверкa aдеквaтности − это оценкa достоверности построенной
мaтемaтической модели, исследование ее соответствия
изучaемому объекту. Для этого вычисляется оценка дисперсии
адекватности:
2
=
Sад
1 N Э
(y j − y Pj )2 ,
∑
N − B j=1
(2.25)
где B − число значимых коэффициентов регрессии. Далее
вычисляют расчетное значение критерия Фишера:
FP =
2
макс{Sад
,S2Y }
2
мин{Sад
,SY2 }
.
(2.26)
По таблице находят табличное значение критерия
Фишера FT [1]. Оно зависит от доверительной вероятности P,
числа степеней свободы fад= N−B и f = N(n−1). На основании
этого делается вывод об адекватности или неадекватности
уравнения регрессии. Уравнение регрессии считается
адекватным, если выполняется условие: Fp ≤ FT .
Чем меньше B, тем больше N-B − в этом одна из главных
целей, достигаемых при исключении незначимых членов.
26
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Если уравнение неадекватно, переходят к более сложной
модели (например, повышают степень многочлена), для чего
обычно требуется постановка добавочных опытов. Иногда
можно обойтись без дополнительного эксперимента, если
соответствующим образом преобразовать переменные у или х.
Интерпретация уравнений регрессии
Интерпретация уравнений регрессии − важнейший этап
моделирования процессов при использовании планирования
эксперимента. Интерпретация включает анализ, прежде всего
влияния отдельных факторов и их взаимодействий, а затем −
особенностей поведения функции отклика в различных частях
изученной области факторного пространства.
Влияние факторов проще всего анализировать по
уравнению 1-й степени. Здесь вначале оценивается знак
коэффициента регрессии, показывающий в какую сторону −
увеличения или уменьшения − влияет на отклик данный фактор.
Таким образом, эмпирический метод составления
математического описания носит итерационный характер,
блочный алгоритм моделирования процесса эмпирическим
методом представлен на рис. 2.4.
Достоинства эмпирического метода построения
математического описания:
• простота описания;
• доступность получения моделей;
• возможность построения модели при отсутствии
теории процесса.
27
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Рис. 2.4. Эмпирический метод
построения математической модели
Недостатки эмпирического метода построения
математического описания:
• отсутствие прогнозируемости в той области, где
эксперимент не проводился; невозможность экстраполяции
результатов;
• необходимость
проведения
большого
числа
экспериментов (часто сложных и дорогостоящих);
• отсутствие ясного физического смысла у параметров
модели;
• непереносимость результатов на другой аппарат.
Главная область применения эмпирического метода
построения
математического
описания
–
системы
автоматизации.
28
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
3. ТЕОРЕТИЧЕСКИЙ МЕТОД ПОСТРОЕНИЯ
МАТЕМАТИЧЕСКОГО ОПИСАНИЯ
Теоретический метод построения математических
моделей процессов химической технологии предполагает
составление
математического
описания
на
основе
фундаментальных законов естественных наук и требует
детального изучения и глубокого понимания физических и
химических закономерностей происходящих явлений.
Важными для построения математического описания
химико-технологических процессов являются следующие
фундаментальные законы:
1. Законы сохранения массы, импульса и энергии,
согласно которым значения этих субстанций при различных
изменениях в изолированной системе должны оставаться
постоянными. Для систем других типов изменение какой-либо
субстанции должно быть равно величине ее внешнего
источника. Математическая запись законов сохранения имеет
форму уравнений балансов, составление которых является
важной частью построения математических моделей.
2. Законы переноса массы, импульса и энергии и законы
химической кинетики, которые определяют плотность потока
этих субстанций. Законы переноса позволяют определить
интенсивность протекающих процессов и, в конечном счете,
эффективность используемых аппаратов.
3. Законы термодинамики, которые формулируют связь
между
тепловыми,
механическими
и
химическими
воздействиями на систему и изменением функций состояния.
Кроме того, законы термодинамики позволяют определить
условия, при которых перенос любой субстанции приходит к
своему завершению (условия равновесия). Знание условий
равновесия необходимо для определения направления
процессов переноса, а также величины движущей силы.
Перечисленные законы составляют теоретическую
основу всех технологических процессов и позволяют получить
их исчерпывающее математическое описание в виде
дифференциальных уравнений с частными производными
второго порядка с соответствующими условиями однозначности
[2]. Это описание позволяет решать как прямую задачу
29
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
поверочного расчета любого аппарата, то есть зная конструкцию
и размеры аппарата, находить поля скорости, давления,
температуры и концентраций в нем, так и обратную задачу
проектного расчета – определять размеры аппарата по
требуемым значениям перечисленных величин на входе и
выходе из него. Проблема заключается лишь в математической
сложности
решения
этих
задач.
Основной
задачей
моделирования является получение на данной теоретической
основе сокращенного математического описания, которое может
быть решено за разумное машинное время на современном
этапе развития математики и средств вычислительной техники.
Далее кратко представим математические выражения для
перечисленных выше законов.
3.1. Теоретические основы
3.1.1. Законы сохранения
Законы сохранения могут записываться применительно
как ко всей системе или ее большим частям (интегральная
форма), так и к элементарным объемам (локальная форма),
которые можно считать точками пространства. В первом случае
получаются алгебраические соотношения балансов, во втором −
дифференциальные уравнения, связывающие изменение
величин, определяющих перенос массы, энергии и импульса в
пространстве. Важным здесь является понятие элементарного
объема, которое имеет смысл только в случае, если существует
изменение характеристик среды в аппарате по какой-либо
пространственной координате. Обычно это некоторый объем,
выделенный в рассматриваемом аппарате, для которого хотя бы
одно из измерений намного меньше соответствующего
масштаба аппарата, а вдоль этого измерения изучаемая
характеристика среды меняется на бесконечно малую величину.
Важно помнить, что математическое описание строится для
расчета
изменения
макрохарактеристик
(температура,
плотность, концентрация и т.д.) рабочей среды во времени и
пространстве, поэтому величина элементарного объема должна
быть намного больше межмолекулярных расстояний. В
противном случае выделенный объем может содержать
30
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
недостаточное
число
молекул
для
того,
чтобы
макрохарактеристики среды в этом объеме не испытывали
значительные флуктуации с изменением числа частиц.
Закон сохранения массы
Масса не может исчезать, либо возникать, то есть
суммарное количество массы в закрытой системе неизменно
(закрытая система не обменивается массой с окружающей
средой), следовательно, ∆М = 0 или dM/dt = 0.
Эти соотношения всегда справедливы для всей массы
системы и в отсутствие химических превращений могут
распространяться на ее компоненты. В противном случае
говорят, что для компонента i существует источник (сток)
массы. Условие закрытой системы не всегда реализуется для
любого выделенного объема внутри реального аппарата, однако
для всего аппарата в целом оно всегда справедливо.
Обобщенное уравнение материального баланса имеет
вид
Приход массы – Расход массы = Накопление массы
Если приход вещества больше расхода, то вещество
накапливается (положительное накопление), если меньше, то
убывает (убыль, или «отрицательное накопление»). В
стационарном режиме обобщенное уравнение запишется в виде
Приход массы = Расход массы.
Следствием закона сохранения массы для элементарного
объема движущейся среды является уравнение сплошности:
∂ρ / ∂t + div ( ρυ) = 0 ,
где ρ – плотность вещества; υ – скорость движения среды.
Закон сохранения энергии
Суть закона сохранения энергии состоит в том, что
энергия не может исчезать, либо возникать, она лишь переходит
из одной формы в другую, то есть dE = 0 или dE/dt = 0, где Е –
31
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
полная энергия системы, которая складывается из внутренней
энергии, кинетической и потенциальной.
Для описания процессов, в которых происходит обмен
энергией в виде тепла, важным является выражение прихода и
расхода тепла через изменение функций состояния системы.
Выбор функции состояния определяется условиями протекания
теплообменного процесса. Например, по второму закону
термодинамики изменение внутренней энергии системы может
происходить за счет теплообмена и совершения работы
dU = δQ − PdV ,
где U – внутренняя энергия системы, которая является частью
полной энергии Е за вычетом энергии движения системы как
целого и ее потенциальной энергии в поле внешних сил. В этом
выражении функциями состояния являются внутренняя энергия,
давление и объем. Количество тепла Q не является функцией
состояния, так как зависит от условий протекания процесса,
поэтому в выражении второго закона термодинамики малое
количество тепла обозначается символом δ . Если процесс
теплообмена протекает при постоянном объеме (V=const), то
количество тепла, переданного системе можно определить как
δQ = dU = C V dT ,
где Cv – изохорная теплоемкость.
В промышленных процессах чаще реализуется ситуация,
когда Р=const. В этом случае за функцию состояния удобно
принять энтальпию, которая имеет следующее определение:
H = U + PV.
Тогда выражение второго закона термодинамики будет иметь
вид
dH = δQ + VdP ,
а количество тепла переданного системе:
δQ = dH = c P dT
или
Q = H1 − H 2 = c P (T1 − T2 ) .
Запишем уравнение теплового баланса для аппарата,
схематично изображенного на рис. 3.1, в котором происходит
теплообмен между теплоносителями через разделяющую их
стенку.
32
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Рис. 3.1. Схема аппарата для составления теплового
баланса
За
время пребывания в аппарате
энтальпии
теплоносителей меняются на величину переданного тепла:
Q = dH1 = dH 2
или
G1c p1 (T1H − T1K ) = G 2 c p2 (T2K − T2H ),
где Gi − массовый расход теплоносителя; Тiн, Тiк – температура
теплоносителя начальная и конечная; срi – теплоемкость
теплоносителя.
Если в движущейся среде существуют внутренние
источники (стоки) тепла, например в результате протекания
химической реакции, то в уравнении теплового баланса они
учитываются соответствующими членами:
G1c p1 (T1H − T1K ) = G 2 c p2 (T2K − T2H ) + Q r
где Q r – количество тепла, выделяемое (поглощаемое) в
результате химической реакции.
Закон сохранения импульса
Суть закона сохранения импульса состоит в том, что
суммарный импульс изолированной системы есть величина
постоянная:
ur
ur
dP
∆P = 0 ,
=0.
dt
Если же система находится под воздействием внешних
сил, то производная от импульса системы по времени равна
результирующей силе, действующей на систему. Следствием
33
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
закона сохранения импульса для элементарного объема
движущейся среды является известная система уравнений
Навье−Стокса.
3.1.2. Законы термодинамики и термодинамического
равновесия
Термодинамика рассматривает явления, обусловленные
поведением большого числа тел (атомов, молекул), из которых
состоят термодинамические системы. В отличие от
молекулярно−статистических теорий в термодинамических
методах изучения равновесных свойств веществ не используют
явно представления о их молекулярном строении. В
термодинамике состояние системы принято описывать в
терминах макроскопических переменных состояния (параметров
состояния), таких как объем V, давление Р, температура Т,
число молей компонентов Nk. Для того что бы задать состояние
системы, достаточно зафиксировать небольшое количество
параметров, все остальные переменные, а также функции
состояния будут при этом иметь строго определенные значения,
поэтому
математически
термодинамические
методы
оказываются намного проще, чем молекулярно-статистические,
где состояние системы определяется 6N переменными (N −
число молекул в системе). Опираясь на феноменологические
законы (начала) в термодинамике получают большое количество
важных для практики следствий, некоторые из которых будут
рассмотрены ниже. Однако необходимо помнить о главном
ограничении термодинамических методов – они позволяют
установить закономерности в поведении свойств веществ, но не
позволяют определить величину самого свойства для данного
вещества.
Формулировка второго начала термодинамики для
открытых систем записывается в виде
n
dU = TdS − PdV + ∑ µi dN i ,
i =1
(3.1)
где S – энтропия; µi и Ni – химический потенциал и число молей
(молекул) компонента i. В соотношении (3.1) внутренняя
34
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
энергия и ее независимые переменные S, V и Ni являются
экстенсивными величинами, то есть они пропорциональны
размерам системы, поэтому их можно представить как: U = Nu ,
S = Ns , V = Nv , N i = Nx i , где N – общее число молей;
величины, обозначенные прописными буквами, являются
удельными, то есть отнесенными к одному молю вещества, xi –
мольная доля i-го компонента. Подставляя эти выражения в (3.1)
и проводя дифференцирование, можно получить важное
следствие:
n
du = Tds − Pdv + ∑ µi dx i
i =1
(3.2)
и
n
n
i =1
i =1
u = Ts − Pv + ∑ µi x i или U = TS − PV + ∑ µi N i .
(3.3)
Если теперь записать дифференциал внутренней энергии,
учитывая (3.3),
n
n
i =1
i =1
dU = TdS + SdT − PdV − VdP + ∑ µi dN i + ∑ N i dµi
и сравнить полученное выражение с (3.1), получим важное и
полезное для практического применения соотношение
Гиббса−Дюгема:
n
SdT − VdP + ∑ N i dµi = 0 ,
i =1
(3.4)
которое показывает, что изменение интенсивных величин
системы (Р, T, µ) не может происходить независимо. При
постоянных Р и T соотношение (3.4) сводится к
n
∑ N dµ
i
i =1
i
=0
или
n
∑ x dµ
i =1
35
i
i
=0.
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Эти выражения часто используются при анализе
фазового равновесия, а также при составлении моделей для
определения химического потенциала. Часто их записывают в
следующем виде:
n
∑x
i =1
i
∂µi
= 0,
∂x k
(3.5)
например, для бинарной смеси получим:
x1
∂µ1
∂µ
+ x2 2 = 0 .
∂x1
∂x1
Соотношение (3.1) является не единственно возможной
математической формулировкой второго закона термодинамики.
Выбор в качестве независимых переменных энтропии не удобен
с практической точки зрения, так как энтропия не может быть
измерена непосредственно в отличие от Р, V, T. Если
термодинамический анализ процесса удобнее проводить в
других переменных, уравнение (3.1) можно переформулировать,
заменив внутреннюю энергию на другую функцию состояния
(термодинамический потенциал), соответствующую выбранным
переменным. Например, для переменных:
V, T, Ni − F = U − TS =
n
∑µ N
i =1
i
i
− PV − свободная энергия;
n
dF = −SdT − PdV + ∑ µi dN i ;
i =1
(3.6)
Р, T, Ni − G = U − TS + PV =
n
∑µ N
i =1
i
i
− энергия Гиббса;
n
dG = SdT + VdP + ∑ µi dN i ;
i =1
(3.7)
Р, S, Ni − H = U + PV = TS +
n
∑µ N
i =1
36
i
i
− энтальпия;
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
n
dH = TdS + VdP + ∑ µi dN i .
i =1
(3.8)
Важным здесь является то, что для получения полной
информации о термодинамических свойствах системы
достаточно знать какой-либо один из термодинамических
потенциалов как функцию своих переменных. Все остальные
термодинамические свойства определяются в этом случае из
разнообразных дифференциальных соотношений Максвелла.
Например, химический потенциал i-го компонента можно
определить как
µi =
∂U
∂N i
=
S,V,N j≠i
∂F
∂N i
=
V,T,N j≠i
∂G
∂N i
=
P,T,N j≠i
∂H
∂Ni
,
S,P,N j≠i
(3.9)
давление получить из соотношения
P=−
∂F
∂U
.
=−
∂V T,Ni
∂V S,Ni
(3.10)
Можно вывести и более сложные дифференциальные
соотношения типа
∂S
∂P
=
и др. [3]. Однако, оставаясь в
∂V T ∂T V
рамках термодинамики, нельзя найти выражения для какоголибо термодинамического потенциала в виде явной функции от
соответствующих переменных, поэтому для определения
термодинамических свойств веществ необходимо привлекать
экспериментальные данные или молекулярно-статистические
методы. Ввиду того что термодинамические потенциалы (U, F,
G, H) нельзя непосредственно экспериментально измерить, на
практике
находят
функцию,
характеризующую
термодинамическое
состояние
системы
между
легко
измеряемыми величинами p, T, V и Ni. Уравнение типа
P = f (T, V, N1 , N 2 ,K N n ) , или уменьшая число переменных на
одну
называют
уравнением
P = f (T, v, x1 , x 2 ,K x n −1 ) ,
состояния. Зная уравнение состояния, и используя
дифференциальные соотношения Максвелла, также можно
37
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
определить
любые
термодинамические
характеристики
системы. Строгую математическую формулировку уравнения
состояния можно получить только в рамках молекулярностатистической теории. Практически активно используются
различные модельные уравнения состояния, некоторые из
которых будут рассмотрены ниже.
Методы термодинамики позволяют сформулировать
условия равновесия многокомпонентных многофазных систем.
Эти условия являются следствием принципа максимума
энтропии в равновесном состоянии системы. При равновесии
гетерогенной системы из Ф фаз с n компонентами температура,
давление и химические потенциалы каждого компонента во всех
фазах одинаковы:
T ' = T '' = T ''' = ... = T ( Ф)
P ' = P '' = P ''' = ... = P ( Ф)
µ1' = µ1'' = µ1''' = ... = µ1( Ф)
µ '2 = µ ''2 = µ '''2 = ... = µ (2Ф)
.....................................
µ 'n = µ ''n = µ '''n = ... = µ (nФ)
(3.11)
Здесь верхний индекс означает фазу, а нижний –
компонент. Данные условия равновесия гетерогенной системы
позволяют определить количество фаз (состоящих из
нескольких компонентов), способных одновременно находиться
в равновесии, или число независимых переменных гетерогенной
системы, которые можно изменять, не нарушая ее равновесия.
Эта задача была решена Гиббсом, поэтому полученный им
результат называется правилом фаз Гиббса: C = n − Ф + 2 , где
С – число степеней свободы, n – число компонентов, Ф – число
фаз системы. Число степеней свободы можно интерпретировать
и как число параметров, которые нужно задать для
однозначного определения состояния системы.
Поскольку химический потенциал играет центральную
роль в описании условий равновесия, а непосредственно
экспериментально его величина не может быть измерена,
38
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
рассмотрим варианты записи химического потенциала и его
определение на основе уравнения состояния.
В однокомпонентном случае химический потенциал
наиболее просто определяется согласно (3.7):
∂G
G
,а
=V.
∂P T,N
N
µ=
Откуда
p
µ
1
∫0 dµ = N ∫0 Vdp .
µ
p
Таким образом, явный вид выражения химического
потенциала получается, если известно уравнение состояния.
Например, для уравнения состояния идеального газа
PV = NRT получим:
p
1
 P 
dP = RT ln  0  .
p
P 
p0
µ(T, P) − µ (T, P ) = RT ∫
0
0
В качестве стандартного состояния (системы отсчета)
удобно взять состояние при P 0 = 1 , тогда µ 0 (T, P 0 = 1) , будем
обозначать просто µ 0 (T) , а химический потенциал запишем как
µ(T, P) = µ 0 (T) + RT ln ( P ) .
(3.12)
В случае смеси химический потенциал i-го компонента
определяется согласно (3.9):
µi (T, P, x) − µi0 (T) =
∂ (G − G 0 )
,
∂N i
P,T,N
j≠i
здесь x обозначает все n-1 концентрации компонентов смеси.
Учитывая, что для смеси уравнение состояния
идеального газа будет PV =
N i RT ,
∑
p
1
 P 
dP = ∑ N i RT ln  0  ,
P
P 
p0
G − G = ∑ Ni RT ∫
0
39
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
а также что за стандартное состояние принимается химический
потенциал чистого компонента, находящегося при единичном
давлении Pi0 = P 0 x i , получим
µi (T, P, x) = µi0 (T) + RT ln ( Pi ) ,
(3.13)
где Pi = Px i − парциальное давление компонента в смеси.
Соотношения (3.12), (3.13) верны только для случая
идеального газа. Если рассматриваемая система не подчиняется
уравнению состояния идеального газа, вводится понятие
фугитивности (летучести) для чистого вещества f = Pγ f или
смеси fi = Pi γ fi . Коэффициент фугитивности γ f зависит от
температуры и давления для чистого вещества, а также еще и от
состава в случае смеси:
µ(T, P) = µ 0 (T) + RT ln ( f ) ;
µi (T, P, x) = µi0 (T) + RT ln ( f i ) .
(3.14)
Очевидно, что для идеального газа γ f = 1 .
Выбор в качестве стандартного состояние идеального
газа, находящегося при единичном давлении, не всегда удобно,
поэтому на практике часто применяются и другие возможные
варианты систем отсчета в выражении химического потенциала
[4]. Например, для компонентов жидкой смеси химический
потенциал обычно определяют как
µi (T, P, x) = µi0 (T, P) + RT ln ( a i ) ,
(3.15)
где µi0 (T, P) − химический потенциал чистого компонента,
находящегося при тех же Т и Р, что и смесь; a i = γ i x i −
активность компонента;
компонента i.
γi
−
коэффициент
активности
fi
, где fi0 − фугитивность чистого
fi0
компонента i при заданных Т и Р. Если γ i = 1 , то такая смесь
Очевидно, что a i =
40
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
называется идеальной. С молекулярно-статистической точки
зрения идеальная смесь – смесь, состоящая из компонентов с
близкими молекулярными характеристиками, то есть размеры
молекул и энергии притяжения, отталкивания у компонентов
смеси близки.
Коэффициент фугитивности может быть определен, если
известно
уравнение
состояния
реального
вещества.
Традиционно уравнение состояния записывают в переменных V
и Т (хотя это не принципиально), поэтому для определения
коэффициента
фугитивности
удобнее
использовать
термодинамическую связь химического потенциала со
свободной энергией. Для чистого вещества согласно (3.6)
получим
ln(Pγ f ) =
(3.16)
µ − µ 0 F − F0 PV − PV 0
=
+
RT
NRT
NRT
Уравнение состояния любого вещества формально
можно представить в виде
PV = ZV NRT ,
где ZV − фактор сжимаемости, функция, которая показывает
отклонение в поведении данного вещества от идеального газа.
Свободная энергия определяется на основе соотношения
(3.10):
V
F − F = − NRT ∫
0
V0
V
 V ZV − 1

ZV
1
dV = − NRT  ∫
dV + ∫ dV 
V
V
V0
 ∞ V

.
Из подынтегрального выражения была явно выделена
идеальная часть, после чего, так как V 0 соответствует
идеальногазовому состоянию, нижний предел интегрирования в
первом члене правой части можно заменить на бесконечность.
Производя интегрирование идеальной части и учитывая,
что
V ZV P 0
=
, окончательно получим
V0
P
41
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
V Z − 1
 Zv  
F − F0 = − NRT  ∫ V
dV + ln 

 P 
∞ V
(3.17)
Подставляя (3.17) в (3.16), получим выражение для
коэффициента фугитивности чистого вещества:
ZV − 1
dV − ln ( ZV ) .
V
∞
V
ln( γ f ) = (ZV − 1) − ∫
(3.18)
Это общее выражение, справедливое для любых веществ.
Коэффициент фугитивности для конкретного вещества
получается
подстановкой
соответствующего
уравнения
состояния. В зависимости от сложности математической
формулировки уравнения состояния интегрирование в (3.18)
выполняется либо аналитически, либо численно.
В многокомпонентном варианте выражение для
фугитивности компонента i определяется следующим образом:
ln(Pi γ fi ) =
µi − µi0
1 ∂ (F − F0 )
=
.
RT
RT ∂Ni V,T,N
j≠i
Зная определение свободной энергии через уравнение
состояния (3.17) несложно записать
V
1 ∂P
ln( γ fi ) = − ∫ 
RT ∂N i
∞

−
V,T,N j≠i
1
 dV − ln ( ZV ) .
V

(3.19)
Если состав смеси представлен в мольных долях, что
обычно удобнее, можно непосредственно дифференцировать по
мольным долям, однако при этом необходима некоторая
осторожность. В отличие от числа молей, мольные доли не
n
могут быть независимыми (
∑x
i =1
i
= 1 ), поэтому при взятии
частной производной в (3.19), по крайней мере две мольные
доли должны изменяться одновременно. Для этого в уравнении
состояния мольную долю одного из компонентов необходимо
выразить через мольные доли остальных компонентов:
42
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
n −1
x k = 1 − ∑ xi .
i≠ k
Кроме того,
n
N i = x i ∑ N j и dNi = dx i N + x i dNi ,
j=1
откуда
dNi =
Ndx i
.
1 − xi
Тогда выражение (3.19) можно переписать в следующем
виде:
1 − x ∂P
i
ln( γ fi ) = − ∫ 
 RTN ∂x i
∞
V
−
V,T
1
 dV − ln ( ZV ) .
V 
Подобным образом, на основе уравнения состояния и
термодинамических соотношений Максвелла можно получить
выражения для любых термодинамических характеристик
рассматриваемой системы.
Приведенные
выше
определения
химического
потенциала и условия термодинамического равновесия фаз дают
возможность записать взаимосвязь равновесных концентраций
компонентов в фазах. Рассмотрим равновесие между паром и
жидкостью. Тогда условие равновесия для i-го компонента
запишется в виде
µi0 (T) + RT ln ( f i ) = µi0 (T, P) + RT ln ( a i )
(3.20)
Или, используя приведенные выше определения,
yi =
f i0 γ i x i
,
γ fi P
где y – концентрация в паре, x − в жидкости.
В инженерных расчетах соотношение (3.20) часто
представляют формально в виде yi = mi x i , однако не следует
забывать, что в общем случае коэффициент распределения m
является сложной функцией параметров состояния и
концентраций компонентов в фазах.
43
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Теперь запишем общую систему уравнений для расчета
фазового равновесия пар−жидкость. По правилу фаз Гиббса для
задания
состояния
двухфазной
системы
необходимо
зафиксировать n (число компонентов) переменных. Обычно
задают состав одной из фаз, например, жидкой, n−1
концентрацией, а также давление или температуру.
Несмотря на эквивалентность этих двух способов
задания состояния системы, трудоемкость расчетного алгоритма
при решении этих задач оказывается неодинаковой. Это
опять−таки связано с традициями представления уравнения
состояния в виде давления как функции от температуры, а не
наоборот, поэтому при задании состояния давлением и составом
требуется либо решение системы большего числа уравнений,
либо использование более сложной итерационной процедуры
для решения этой системы.
Пусть известны n−1 концентраций компонентов в
жидкой фазе и температура T, одинаковая во всех фазах. Тогда
из условий термодинамического равновесия следует следующая
система уравнений:


P(T, vп , y1 , y 2 ,K y n −1 ) = P(T, v ж , x1 , x 2 ,K x n −1 )

fi0 γ i x i
y
=
,
i=1K n
 i
γ
P
fi

n
∑ yi = 1
 i=1
(3.21)
Решение системы n+2 уравнений позволяет определить n
концентраций в паровой фазе, удельные объемы пара и
жидкости, а также независимо по уравнению состояния
давление в фазах.
Если исходными данными являются n−1 концентраций
компонентов в жидкой фазе и давление Р0, то уравнение
состояния для определения давления необходимо включить в
общую систему уравнений:
44
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
P(T, vп , y1 , y 2 ,K y n −1 ) = P0
P(T, v , x , x ,K x ) = P
ж
1
2
n −1
0

0

fi γ i x i
,
i=1K n
 yi =
γ
P
f

i
n
∑ yi = 1
 i=1
(3.22)
В итоге получим систему n+3 уравнений, решение
которой дает значения n концентраций в паровой фазе,
удельные объемы пара и жидкости и температуру.
Решение задачи расчета условий равновесия фаз на
основе уравнения состояния чрезвычайно выгодно, так как
помимо равновесных составов фаз, позволяет получить полную
информацию о термодинамических свойствах системы.
Например, по результатам решения систем (3.21) или (3.22) и
уравнению состояния не составляет большой сложности
определить многие важные для анализа тепломассообменных
процессов свойства: плотность пара и жидкости, энтальпию фаз
и теплоту парообразования и т.д. Однако практическая
реализация расчета условий фазового равновесий на основе
уравнения состояния используется достаточно редко. Основной
проблемой здесь является отсутствие надежных уравнений
состояния, способных с удовлетворительной точностью
описывать как паровое, так и жидкое состояние
многокомпонентных
неидеальных
систем,
поэтому
в
инженерных
методах
отдельные
задачи
определения
термодинамических свойств многофазных многокомпонентных
систем
разделяют
на
задачи
расчета
различных
термодинамических характеристик. Для этих характеристик
приходится создавать различные модели эмпирического или
полуэмпирического типа, которые имеют ограничения в
применении, как по области термодинамических состояний, так
и по веществам и их смесям.
Задачу расчета характеристик фазового равновесия
пар−жидкость, можно разделить на следующие части:
определение равновесного состава фаз и давления или
45
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
температуры, определения плотности фаз и определение
различных важных для анализа теплофизических свойств
(теплоемкость, энтальпия, теплота парообразовании и т.д.). Для
определения состава фаз и давления достаточно следующей
системы уравнений:

fi0 γ i x i
=
y
,
 i
γ fi P

n
 y =1
i
∑
i =1
i=1K n
(3.23)
Решение этой системы возможно, если дополнить ее
моделями определения коэффициентов фугитивности и
активности в зависимости от Т, Р и состава. Часто систему
уравнений (3.23) можно упростить, если принять, что паровая
фаза подчиняется уравнению состояния идеального газа. Тогда
γ fi = 1 и при небольших давлениях значение фугитивности
близко к давлению насыщенных паров чистого компонента
fi0 (T, p) ≈ Pi0 (T) [5]. В итоге получим

Pi0 γ i x i
=
y
,
 i
P
n
 y =1
i
∑
i =1
i=1K n
(3.24)
Это наиболее часто используемая формулировка расчета
состава равновесных пара и жидкости. Систему уравнений
(3.24) требуется дополнить зависимостью коэффициента
активности от Т, Р и состава и зависимостью давления
насыщенных паров чистых компонентов от температуры.
Некоторые модели коэффициентов активности (Вильсона,
НРТЛ, ЮНИФАК) будут рассмотрены далее, а подходящий вид
аппроксимации для давления насыщенных паров получают из
известного
термодинамического
соотношения
Клаузиуса−Клайперона:
46
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
∂P 0 (T) ∆H
=
,
∂T
T∆V
(3.25)
здесь ∆H = H П − H Ж , ∆V = VП − VЖ − разность энтальпий и
объемов паровой и жидкой фаз. Используя формальный вид
уравнения состояния, перепишем (3.25) в виде
∂ ln ( P 0 (T) )
1
∂ 
T
=
∆H
,
R∆ZV
(3.26)
где ∆ZV − разность факторов сжимаемости паровой и жидкой
фаз.
Далее экспериментальные данные показывают, что для
многих веществ, отношение
∆H
можно аппроксимировать
∆ZV
несложными зависимостями. Например, если принять это
отношение постоянной величиной, то интегрирование (3.26)
приводит к уравнению Клайперона:
ln(P 0 ) = A −
B
,
T
(3.27)
где А и В − параметры модели, которые определяются по
имеющимся экспериментальным данным, например методом
наименьших квадратов.
Также можно предложить модели для давления
насыщенных паров от температуры с большим числом
параметров:
ln(P 0 ) = A −
B
− уравнение Антуана;
T+С
ln(P 0 ) = A −
B
+ DT 6 + C ln(T) − уравнение Риделя.
T
(3.28)
(3.29)
47
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Таким образом, задача расчета равновесных составов фаз
оказывается замкнутой. Если известен состав жидкой фазы и
температура, то алгоритм решения системы уравнений (3.24)
следующий:
1. Для заданной температуры рассчитывают давления
паров чистых компонентов Pi0 (T) , например по какому-либо из
уравнений (3.27) − (3.29).
2. Для заданного состава и температуры по выбранной
модели рассчитывают коэффициенты активности компонентов
γi. Несмотря на то, что для определения коэффициента
активности необходимо задать состав, температуру и давление,
на практически часто используемые модели (Вильсона, НРТЛ,
ЮНИФАК), зависимость от давления не учитывают,
предполагая ее малой, поэтому в хорошем приближении γi
может быть определена как функция только состава и
температуры.
3. По (3.24), рассчитывают парциальные давления
компонентов в смеси Pi = Pi0 (T) γ i x i и определяют полное
n
давление
∑P = P .
i =1
i
4. Определяют равновесные концентрации в паровой
фазе: yi =
Pi
.
P
Если исходными данными для расчета являются состав
жидкой фазы и давление P0, то в силу нелинейности уравнений
алгоритм решения системы уравнений (3.24) будет носить
итерационный характер:
1. Задают начальное приближение температуры Т0.
2. Производятся расчеты с п.1 по п.3 согласно
предыдущему алгоритму.
3. Проверяют выполнение условия P0 − P < ε , где ε
некоторое малое число, определяющее точность расчета. Если
это условие выполняется, то переходят к следующему пункту, в
противном случае задают новое приближение температуры Т0 и
расчет повторяют с п. 2.
48
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
4. Определяют равновесные концентрации в паровой
фазе: yi =
Pi
.
P
Если для анализа процесса требуется знание плотностей
и других термодинамических характеристик равновесных фаз,
то для их расчета необходимо дополнительно привлекать
соответствующие модели.
Расчетные схемы фазовых равновесий других типов
(газ−жидкость, жидкость−жидкость) строятся по аналогичной
схеме.
Для фазовых равновесий типа пар(газ)−жидкость
уравнения, связывающие концентрации в фазах, могут быть еще
более упрощены. Например, если жидкая смесь близка к
Pi0 x i
, что является законом Рауля.
идеальной ( γ i = 1 ), то yi =
P
Закон Рауля справедлив и для неидеальных смесей, но только
для компонентов, концентрации которых x i → 1 . Другой
асимптотический случай x i → 0 , обычно имеющий место в
системах газ−жидкость (газы плохо растворяются в жидкостях),
хорошо описывается законом Генри
yi =
Ex i
, где Е –
P
коэффициента Генри.
Химическое равновесие
Химическое равновесие – это термодинамическое
равновесие в системе, между компонентами которой происходят
химические реакции.
Рассмотрим химическую реакцию, происходящую в
закрытой системе
νA A + νBB
ν CC ,
где νА, νВ, νС – стехиометрические коэффициенты участвующих в
реакции компонентов А,В,С соответственно.
Для анализа химических превращений в термодинамике
вводится новая переменная состояния
– химическое сродство
[6].
49
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
= (ν A µ A + ν Bµ B − ν C µ C ) .
В состоянии равновесия химическое сродство реакции
равно нулю, то есть
= 0 . Таким образом, можно заключить,
что в состоянии термодинамического равновесия химические
потенциалы веществ А, В и С принимают такие значения, для
которых
= (ν A µ A + ν Bµ B − ν C µ C ) = 0 .
Используя активность, получим:
ν A (µ 0A (T) + RT ln(a A )) + ν B (µ0B (T) + RT ln(a B )) −
−ν C (µ0C (T) + RT ln(a C )) = 0
0
0
0
ν A µ A +ν Bµ B −ν C µC
a CνC
RT
=
= K(T)
e
νA νB
aA aB
(3.30)
В состоянии химического равновесия система
характеризуется константой равновесия K(T), определяемой
выражением (3.30). Как видно из этого уравнения K(T) −
функция только температуры Т.
Помимо термодинамического определения химического
равновесия
существует
и
кинетическое
определение:
химическое равновесие устанавливается при условиях, что
скорость прямой реакции (исходные реагенты → продукты) r1
становится равной скорости обратного превращения (продукты
→ исходные реагенты) r2:
r1 = r2 .
Вспомним, что под скоростью гомогенной химической
реакции
понимают
изменение
количества
вещества,
вступающего в реакцию или образующегося при реакции за
единицу времени в единице объема системы:
r=
∆C
,
∆τ
где r − скорость реакции; ∆C − изменение концентрации; ∆τ −
интервал времени, в течение которого это изменение
произошло.
Скорости реакции при теоретическом анализе
предпочтительнее выражать через активности, а не
50
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
концентрации, так как состояние равновесия непосредственно
связано с активностями, а не с концентрациями:
r1 = k1a νAA a νBB
r2 = k 2 a CνC
(3.31)
где k1, k2 − константы скорости прямой и обратной реакций.
Сравнивая выражение (3.31) и константу равновесия
(3.30), видим, что
K (T ) =
k1
k2 .
Таким образом, константа равновесия может быть также
связана с константами скоростей, когда скорости реакций
выражены через активности.
3.1.3. Законы переноса и химической кинетики
Когда макросистема находится в равновесии, все ее
термодинамические параметры постоянны по всему объему
системы. Если систему вывести из равновесия и предоставить
самой себе, то она постепенно вернется в равновесное
состояние. При этом будут протекать необратимые процессы,
называемые процессами переноса.
Процессы переноса − это необратимые процессы
пространственного переноса массы, импульса, энергии.
Причины этих процессов – пространственные неоднородности
полей концентраций, скоростей и температур. Перенос
происходит в направлении, обратном градиенту концентрации,
скорости, температуры, что приближает систему к равновесию.
Процессы переноса в покоящейся среде осуществляются
только в результате хаотичного движения молекул
(молекулярный перенос). В текущих средах к этому механизму
переноса добавляется конвективный перенос, а при высоких
числах Рейнольдса еще и турбулентный перенос, связанный с
хаотичным перемещением турбулентных вихрей. Общая теория
процессов переноса изложена в [7], здесь приведем сводку
выражений, описывающих перенос субстанции за счет
молекулярного, конвективного и турбулентного механизмов в
51
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
случае одномерного движения потока вдоль оси Х со скоростью
υx .
Молекулярный механизм:
для потока массы
ur
ur
ji = − Dij ∇Ci
для потока тепла
r
ur
q = −λ∇T
для потока импульса
τ xy = −µ
dυx
dy
Конвективный механизм:
для потока массы
r
r
ji = Ci υx
для потока тепла
r
r
q = ρc p T υx
для
импульса
потока
τ xx = ρυx υx
Турбулентный механизм:
для потока массы
r
ur
ji = − DT ∇Ci
для потока тепла
r
ur
q = −λ T ∇T
для потока импульса
τ xy = −µ T
dυx
dy
где j, q, τ – поток массы, тепла и импульса; υ − конвективная
скорость; Dij – коэффициент бинарной диффузии; Сi –
концентрация компонента i; ρ − плотность, λ − коэффициент
молекулярной
теплопроводности;
µ
−
коэффициент
динамической молекулярной вязкости; DТ – коэффициент
турбулентной диффузии; λТ − коэффициент турбулентной
теплопроводности;
µТ
−
коэффициент
динамической
турбулентной вязкости.
Заметим, что записанные выражения для потока массы
молекулярным и турбулентным механизмом будут строго
справедливы только для двухкомпонентной системы [7].
Нетрудно убедиться в аналогии этих уравнений.
Конвективный поток представляет собой произведение
r
переносимой субстанции в единичном объеме (ρ, ρЕ, ρ υ ) на
конвективную скорость. Потоки за счет молекулярного или
турбулентного
механизмов
есть
произведение
соответствующего коэффициента переноса (D, λ, µ) на
52
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
движущую силу процесса. Следует обратить внимание, что если
D = f (T, P, x) , то DT = f (υ, T, P, x) .
При определении скорости химической реакции
различают гомогенный (в объеме) и гетерогенный (на
поверхности) варианты ее протекания. Тогда под скоростью
химической реакции понимают изменение массы вещества,
происходящее в единицу времени в единице объема для
гомогенной и на единицу поверхности для гетерогенной
реакции. Для реакции типа
ν A A + ν BB → ν CC
в гомогенном варианте:
rA =
1 dN A
V dτ
или, если объем системы остается постоянным,
rA =
dC A
.
dτ
В гетерогенном варианте:
rA =
1 dN A
.
F dτ
Подобным образом можно определить скорость
химической реакции через изменение массы как компонента В,
так и компонента С. Очевидно, что в общем случае rA ≠ rB ≠ rC ,
поэтому удобнее скорость химической реакции определить как
r=
rA
=
νA
rB
νB
=
rC
νC
.
Так же вводят понятие степени полноты реакции,
которое равно отношению изменения числа молей N 0i − Ni
любого из участвующих в реакции веществ к его
стехиометрическому коэффициенту v:
ξ=
N 0A − N A
νA
=
N 0B − N B
νB
=
N 0C − N C
νC
.
где N 0i , N i − начальное и конечное число молей i − го вещества.
Тогда:
53
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
r=
1 dξ
V dτ
r=
1 dξ
.
F dτ
или
Выражение для нахождения скорости химической
реакции можно получить из простых рассуждений. Химическое
превращение происходит при столкновении молекул реагентов.
Число столкновений зависит от концентрации компонентов. В
тоже время не каждое столкновение приводит к химическому
процессу, а лишь та их часть, которая происходит при
определенных кинетических энергиях молекул. Тогда
формально скорость химической реакции можно принять
пропорциональной концентрациям реагентов:
r = kCnAA CBn B ,
(3.32)
где nA, nB – порядок реакции по компоненту А и В
соответственно.
Коэффициент пропорциональности k – константа
скорости химической реакции − зависит от температуры. В
большинстве случаев эта зависимость хорошо описывается
уравнением Аррениуса:
k = A exp(− E / RT) ,
(3.33)
где А – параметр, учитывающий частоту столкновений; Е –
энергия активации.
Уравнения подобного типа (3.32) называют уравнениями
формальной химической кинетики. По своей сути это некоторая
математическая модель, в которой константа скорости и
показатели степени у концентраций являются параметрами,
определяемыми из экспериментальных данных.
Сумма показателей степени называется порядком
реакции. Для простых реакций показатель совпадает со
стехиометрическим
коэффициентом,
для
сложных,
многостадийных реакций такое совпадение не обязательно, а
вид самого уравнения может быть существенно сложнее.
54
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
3.2. Исчерпывающее описание процессов химической
технологии и типовые модели структуры потоков
Основной задачей моделирования процессов химической
технологии является изучение поведения характеристик рабочих
сред во времени и пространстве внутри аппарата. К наиболее
важным характеристикам относятся: температура, давление,
вектор скорости, концентрации компонентов в фазах.
Изменения этих характеристик будут зависеть от их значений на
входе в аппарат, внешних факторов (теплообмен через стенку
аппарата, ввод дополнительных потоков), а также физических и
химических процессов, протекающих внутри аппарата.
Особенно важно, чтобы математическая модель правильно
описывала поведение характеристик среды на выходе из
аппарата, поэтому в некоторые математические модели
(например, ячеечная модель) заведомо закладывается
качественно нереализуемое поведение характеристик рабочей
среды в пространстве аппарата. Тем не менее, на выходе из
аппарата
значение
исследуемой
величины
должно
соответствовать экспериментальным результатам.
Теоретический, или структурный, метод построения
математической модели заключается в составлении уравнений
балансов для элементарного объема в движущейся однородной
среде.
Величина выделенного объема может быть различной,
минимально он должен быть таким, что среда, заключенная
внутри него, еще может считаться сплошной и к ней применимо
понятие термодинамической системы. Максимально размеры
его границ могут достигать масштаба аппарата, если вдоль
направления этой грани не происходит каких−либо заметных
изменений характеристик.
Для примера рассмотрим элементарный объем dV =
dxdydz в движущейся среде, внутри которого имеется источник
химической реакции и происходит межфазный перенос массы
j12. Перенос тепла отсутствует. Для упрощения задачи примем,
что массовый поток j изменяется только вдоль оси Х (рис. 3.2), а
межфазный
перенос
осуществляется
через
нижнюю
поверхность.
55
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Рис. 3.2. Изменение массового потока, направленного
вдоль оси Х и проходящего через элементарный объем dV
Запишем уравнение материального баланса для этого
элементарного объема dV:
jdS + j12 dF − ( j + dj)dS − ri dV =
∂Ci
dV ,
∂τ
(3.34)
где ri = kCin − скорость химического превращения; k –
константа скорости реакции; n – порядок химической реакции.
Поделим уравнение (3.34) на dV, получим
(3.35)
dF
∂Ci ∂ji
+
= j12
− ri
∂τ ∂x
dV
В уравнении (3.34) поток массы ji осуществляется за счет
конвективного и молекулярного механизмов переноса вещества,
который определяется как
ji = υCi − D
dCi
,
dx
а межфазный поток массы j12, можно определить уравнением
массопередачи
j12 = K(Ci − C*i ) ,
где К – коэффициент массопередачи; C*i
концентрация в фазе.
56
− равновесная
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Принимая во внимание, что
dF
= a − это удельная
dV
поверхность контакта фаз, запишем уравнение (3.34) для
одномерного профиля концентраций, учитывающее межфазный
перенос массы и химическое превращение:
(3.36)
∂Ci
∂C
∂ 2C
+ υ i = D 2i + Ka(Ci − C*i ) − kCin .
∂τ
∂x
∂x
В уравнении (3.36) ii = Ka(Ci − C*i ) − kCin обычно
называют источником (стоком) массы.
В реальных аппаратах движение всегда носит
трехмерный характер, поэтому данное уравнение необходимо
дополнить производными по координатам Y и Z и
соответствующими проекциями вектора скорости. В этом
случае для определения проекций вектора скорости на оси
координат нужно решать систему уравнений Навье−Стокса и
уравнения неразрывности.
Подобное описание, представляющее собой систему
дифференциальных уравнений с частными производными
второго порядка, является исчерпывающим для моделирования
процессов переноса массы. К сожалению, получить строгое
аналитическое решение этой системы уравнений нельзя.
Возможности вычислительной техники для численного решения
подобных систем дифференциальных уравнений в настоящее
время оказываются во многих случаях недостаточными,
поэтому на практике применяются модели, имеющие более
простой математический вид. Получить эти модели для
конкретной ситуации можно из исчерпывающего описания
посредством его сокращения на основе оценки масштабов его
членов, например, так как это делается в теории пограничного
слоя. Также математическую модель можно получить,
рассматривая данное явление на основе намного упрощенной
физической картины. Например, для колонных аппаратов,
высота которых больше диаметра ( L D ), можно
предположить, что неравномерности характеристик потока в
поперечном направлении незначительны. Тогда для расчета
профиля концентраций по высоте колонны будем использовать
57
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
уравнение (3.36). При этом влияние переноса субстанции в
других направлениях можно учесть введением некоторых
параметров, значения которых зависят от конструкции аппарата
и условий проведения процесса. Значения этих параметров
можно определить на основе сопоставления результатов расчета
с данными физического эксперимента.
Таким образом, получается диффузионная модель,
согласно которой структура потока описывается уравнением,
аналогичным уравнению (3.36), но вместо коэффициента
молекулярной диффузии вводится эффективный коэффициент
продольного
перемешивания
DL.
В
результате
однопараметрическая диффузионная модель будет иметь вид
(3.37)
∂Ci
∂C
∂ 2 Ci
+ υ i = DL
+i .
∂τ
∂x
∂x 2
Уравнение (3.37) описывает нестационарный режим.
Стационарный режим характеризуется установившимися
значениями параметров, поэтому для его описания слагаемое,
учитывающее изменение концентрации вещества во времени,
приравнивается к нулю:
∂Ci
=0.
∂τ
В правой части уравнения (3.37) первое слагаемое
учитывает перенос вещества за счет некоторого эффективного
диффузионного потока вдоль направления движения среды.
Для получения однозначного решения уравнения (3.37)
должны быть заданы одно начальное и два граничных условия.
В качестве начального условия задается профиль концентраций
по аппарату в начальный момент времени:
C0i = Ci (τ0 , x) .
(3.38)
Граничные условия обычно записываются в неявном
виде, предложенном Данквертсом (рис. 3.3):
рассматривая условие материального баланса на входе в
аппарат, получим:
58
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
при x = 0
υ(Ci − Cвх i ) = − DL
(3.39)
dCi
;
dx
на выходе из аппарата принимают
при x = L dCi / dx = 0 .
(3.40)
Рис. 3.3. К постановке граничных условий
Часто
безразмерной
переменные
диффузионную
модель
представляют
в
форме. Для этого вводятся безразмерные
X=x/L; τ =
Vап SL L
=
= ; θ =τ / τ .
G Sυ υ
Подставив их в уравнение (3.37), получим
1 ∂Ci υ ∂Ci D L ∂ 2 Ci
+
= 2
+i
2
τ ∂θ L ∂X L ∂x
или
1 ∂ 2 Ci
∂Ci ∂Ci
+
=
+ iτ ,
∂θ ∂X Pe L ∂X 2
(3.41)
где множитель Pe L = (υL) / D L в левой части уравнения (3.41)
представляет собой критерий Пекле для обратного
перемешивания.
Наряду
с
рассмотренной
однопараметрической
диффузионной моделью используется двухпараметрическая
диффузионная модель. Отличие ее состоит в том, что
перемешивание потока учитывается как в продольном, так и в
радиальном направлении.
59
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Для описания радиального переноса используется
эффективный
коэффициент
диффузии
в
радиальном
направлении в цилиндрической системе координат DR.
Остальные обозначения аналогичны предыдущей модели:
(3.42)
∂Ci
∂C
∂ 2 Ci
∂ 2 Ci 1 ∂Ci
+ υ i = DL
+
D
(
+
)+i .
R
∂τ
∂x
∂x 2
∂r 2 r ∂r
К граничным условиям однопараметрической
диффузионной модели добавляются условия при r = 0 и r = R.
при x = 0 υ(Сi − C i вх ) = DL
dCi
;
dx
dCi
= 0;
dx
∂Сi
= 0;
при r = 0
∂r
∂Сi
при r = R
= 0.
∂r
при x = L
Идя по пути упрощения физической картины явления,
можно предположить отсутствие продольного перемешивания,
тогда из уравнения (3.37) получаем модель идеального
вытеснения (МИВ), согласно которой все элементы потока
движутся по параллельным траекториям с одинаковыми
скоростями. Время пребывания в аппарате для всех элементов
такого потока одинаково.
Математическое выражение для МИВ получим,
приравняв коэффициент продольного перемешивания к нулю в
уравнении (3.37)
(3.43)
∂CA
∂C
+υ A =i.
∂τ
∂x
Начальные и граничные условия будут иметь вид:
начальные условия: τ = 0; C = C(x) .
граничные условия: x = 0; C = Cвх .
МИВ означает идеальный случай, когда любое
перемешивание отсутствует. Другой крайне идеальной моделью
60
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
является модель идеального смешения (МИС), согласно которой
концентрация одинакова по всему объему аппарата и равна
соответствующему значению на выходе. Тогда производная
концентрации по координате в уравнении (3.43) заменяется
разностью концентраций на входе и выходе, отнесенной к
среднему времени пребывания τ :
(3.44)
∂CA (CA − C A вх )
+
=i.
∂τ
τ
Для стационарного режима МИС примет вид
(CA − C A вх )
τ
(3.45)
= i.
Таким образом, статическая МИС представляет собой
алгебраическое уравнение.
Структуру потоков в реальном аппарате можно описать с
помощью ячеечной модели, в рамках которой реальный аппарат
мысленно разбивается на некоторое число одинаковых
последовательно соединенных ячеек (аппаратов) идеального
смешения Vап/m. Количество ячеек m является параметром,
характеризующим ячеечную модель реального аппарата. Если m
= 1, то ячеечная модель переходит в модель идеального
смешения, а в случае m=∝
∝ − в модель идеального вытеснения.
В ряде случаев в аппаратах действительно можно
выделить участки по ходу потока, в каждом из которых режим
близок к идеальному смешению − каскад реакторов идеального
смешения, тарельчатая барботажная колонна, реактор с
кипящим слоем.
При составлении ячеечной модели используются
следующие допущения:
- в каждой ячейке поток имеет структуру идеального
перемешивания и концентрация C не изменяется в пределах
соответствующей ячейки (j = 1, 2, … n − индекс ячейки);
- перемешивание между ячейками отсутствует;
•
- объемный расход V не изменяется.
61
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Кроме того, чаще всего для удобства расчетов
принимается, что:
- объемы каждой из предполагаемых m ячеек одинаковы
и равны Vj;
- сумма объемов всех ячеек равна общему объему зоны,
для которой справедлива ячеечная модель (Vaп = mVj);
- среднее время пребывания частиц в каждой ячейке −
•
τ j = τ / m , а среднее время пребывания в системе − τ = Vап / V .
Поскольку в каждой ячейке соблюдается режим
идеального смешения, для любой j−ой ячейки справедливо
уравнение МИС:
∂C j
∂τ
+
(C j−1 − C j )
τj
= ij ,
(3.46)
где j = 1, 2, ..., m.
Система уравнений представляет собой материальные
балансы в каждой из принятых m ячеек.
Начальные условия и ограничения записываются для
каждой ячейки аналогично модели МИС.
3.2.1. Импульсный ввод индикатора для определения
параметров типовых и комбинированных моделей
структуры потоков
Рассмотренные выше модели описания процессов
массопереноса и химического превращения основаны на
некоторой упрощенной картине структуры потоков в аппарате.
Если принимается не идеальная структура потоков (МИС,
МИВ), то модель обязательно будет иметь один или несколько
настраиваемых параметров (диффузионная – число Pe, ячеечная
– число ячеек m). Значения этих параметров можно определить
только из экспериментальных данных, при этом методы
извлечения значений параметров и необходимые для этого
экспериментальные данные могут быть различны.
Пусть, например, требуется определить параметр
диффузионной модели для описания профиля температур
одного из теплоносителей в теплообменном аппарате. Для
62
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
решения этой задачи можно воспользоваться методом
наименьших квадратов. В этом случае необходимо
экспериментально измерить необходимые характеристики
теплоносителя в нескольких точках по длине аппарата и далее
определить значение параметра, минимизируя разницу
расчетных и экспериментальных данных.
Этот метод в большинстве случаев оказывается
неудобным, так как:
• затруднительно определить изменения характеристик
процесса по длине аппарата;
• необходимо многократно находить решение по
модели, что может быть трудоемким, в случае если модель
имеет вид дифференциального уравнения.
На практике для определения параметров основных
моделей структуры потока широкое распространение получил
метод импульсного ввода индикатора (иногда его называют
методом меченых объемов). Суть этого метода состоит в том,
что в основной поток на входе в аппарат вводится индикатор,
элементарные объемы которого должны двигаться аналогично
элементам основного потока. Фиксируя изменение во времени
концентрации индикатора на выходе из аппарата, получают
кривую отклика, которая может быть приведена в соответствие
с функцией распределения времени пребывания элементов
потока в аппарате. Этот метод лишен перечисленных выше
недостатков.
Вид кривой отклика на импульсное возмущение уже
качественно позволяет судить о близости структуры потоков в
данном аппарате к идеальной.
Например, для предельных случаев, когда отсутствует
любое перемешивание (МИВ), и для идеального перемешивания
(МИС), кривые отклика будут иметь вид, представленных на
рис. 3.4.
На рис. 3.5 представлен график зависимости
концентрации С от времени τ для реального аппарата.
63
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
а
б
Рис. 3.4. Зависимость концентрации С от времени τ: а –
для МИВ; б – для МИС
Рис. 3.5. Зависимость концентрации С от времени τ для
реального аппарата
Таким образом, чем уже пик на кривой отклика, тем
меньше перемешивание, и наоборот.
Количественно значения параметров модели могут быть
определены на основе статистического анализа кривых отклика.
Несложно показать, что кривая отклика связана с функцией
распределения по времени пребывания частиц в аппарате:
64
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
f(τ) =
dN(τ)
=
Ndτ
C(τ)
∞
,
(3.47)
∫ C(τ)dτ
0
где dN(τ) – количество элементов потока, время пребывания
которых в аппарате составляет от τ до τ+dτ; N – общее число
выделенных элементарных объемов в аппарате; М – масса
введенного индикатора; С(τ) – концентрация индикатора на
выходе из аппарата. Условие нормировки:
∞
∫ f(τ)dτ =1 ,
(3.48)
0
то есть вероятность того, что частица, вошедшая в аппарат,
когда-нибудь выйдет из него, равна единице.
Расчет распределения времени пребывания частиц
потока основан на статистическом понятии моментов и связан с
распределением плотности вероятностей. Основные свойства
распределения
случайной
величины
можно
описать
несколькими
числовыми
характеристиками,
которые
определяют
наиболее
существенные
особенности
распределения. Такой системой характеристик являются
моменты распределения случайной величины, которые
систематизируются по трем признакам: по порядку n момента;
по началу отсчета случайной величины; по виду случайной
величины.
Порядок момента может быть любым целым числом.
Практически же рассматривают моменты нулевого, первого,
второго, третьего и четвертого порядков, то есть n = 0, 1, 2, 3, 4.
В зависимости от начала отсчета случайной величины
различают начальные и центральные моменты. Общий вид
начальных моментов фукции распределения таков:
∞
M n = ∫ τ n f (τ)dτ .
0
Каждый из моментов имеет определенный физический
смысл. Нулевой момент – это площадь под кривой; первый
момент характеризует среднее значение (среднее время
пребывания), или математическое ожидание случайной
величины времени пребывания.
65
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Случайные
величины,
отсчитываемые
от
математического ожидания, называются центрированными.
Моменты
центрированной
величины
называются
центральными. Общий вид центральных моментов таков:
∞
µ n = ∫ (τ − τ) n f (τ)dτ .
0
Второй центральный момент характеризует рассеяние
случайной величины относительно среднего времени
пребывания, он называется дисперсией и обозначается через
σ2τ :
∞
σ2τ = µ 2 = ∫ (τ − τ)2 f (τ)dτ .
0
Из
соображения
удобства
часто
используют
безразмерные время пребывания, функцию распределения и
безразмерные моменты:
(3.49)
θ=τ/ τ .
f * (θ) = τ f(τ) .
σθ2 =
σ
τ
2
τ
2
(3.50)
.
Оказывается, параметры многих моделей, включая
диффузионную, ячеечную и их комбинации, можно связать с
моментами. Процедура определения такой связи достаточно
трудоемка [8] и заключается в аналитическом интегрировании
нестационарного решения для данной модели по времени. Для
диффузионной и ячеечной модели можно получить следующие
соотношения, связывающие дисперсию времени пребывания
элементов потока в аппарате σθ2 с параметрами этих моделей.
Ячеечная модель:
Диффузионная модель:
2 2 (1 − e
σθ =
−
Pe
Pe2
− Pe
1
= σ θ2 .
m
2
).
Таким образом, поиск параметров основных моделей
структуры потока методом импульсного ввода индикатора
можно представить следующим алгоритмом:
66
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
• проводится эксперимент методом импульсного ввода
индикатора;
• определяется кривая отклика C( τ) ;
• по кривой отклика находится функция распределения
и необходимые моменты;
• находятся параметры модели.
Ячеечная модель с обратными токами
Ячеечная модель не всегда обеспечивает адекватное
воспроизведение структуры потоков в реальном аппарате. В
связи с этим разработаны модификации такой модели. Одной из
наиболее распространенных модификаций является ячеечная
модель с обратными потоками (рис. 3.6). Согласно этой модели
аппарат рассматривают как последовательность зон с
сосредоточенными параметрами, причем каждая из зон
эквивалентна
ячейке
идеального
перемешивания.
Модифицированная ячеечная модель (рециркуляционная) может
учитывать сложности, возникающие в различных аппаратах при
оценке существования обратных потоков в аппарате. Меняя
параметры данной модели (N – число ячеек, f – доля обратного
потока), можно описывать гидродинамические режимы, близкие
к режимам идеального смешения и вытеснения, а также
режимы,
соответствующие
диффузионной
однопараметрической и ячеечной моделям.
Рис. 3.6. Схема потоков по ячеечной модели с обратными
потоками: L – поток вещества; е – обратный поток вещества; Сj
– концентрация на выходе i-ой ячейки
Запишем материальный баланс для каждой из ячеек с
учетом обратных потоков между ними, если величины обратных
потоков из каждой ячейки равны
67
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
− первая ячейка:
V1
•
dC1 •
= V C Н + eC2 − (V + e)C1 ,
dτ
•
где V1 − объем ячейки, V − объемный расход; е – объемный
расход обратного потока.
− i-я ячейка:
Vj
dC j
dτ
•
•
= (V + e)C j−1 + eC j+1 − (V + e)C j − eC j .
Введем следующие величины: f =
e
•
− доля обратного
V
потока; x =
f
1+ f
Значения параметров m и f определяются в результате
совместного решения следующих уравнений:
M θ2 = 1 +
m(1 − x 2 ) − 2x(1 − x m )
,
m 2 (1 − x) 2
(3.51)
где M θ − безразмерный центральный момент;
M 3θ = 1 +
2 6x(1 + 3x m ) + 3m(1 − x 2 ) 12x(1 + x)(1 − x m )
+
−
.
m2
m 2 (1 − x) 2
m 3 (1 − x 3 )
(3.52)
∫
В этих уравнениях M θm = θm f (θ)dθ .
Комбинированные модели
Не все реальные процессы удается описать при помощи
рассмотренных ранее моделей. В частности, это относится к
процессам, имеющим циркуляционные и байпасные потоки
либо застойные зоны. В таких случаях используются
комбинированные модели.
Особенно сложные потоки часто удобно описывать
комбинированными моделями, построенными как совокупность
зон, соединенных последовательно, параллельно или в
68
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
каком−либо более сложном порядке, структуры потоков
которых различны: зона поршневого потока (идеального
вытеснения); зона потока с идеальным перемешиванием; зона с
продольным перемешиванием; застойная зона. Помимо этого,
могут
наблюдаться
локальные
потоки:
байпасный,
циркуляционный, проскальзывания и т.д.
Следует иметь в виду, что увеличением числа зон можно
описать процесс любой сложности, но математическое
моделирование при этом усложняется.
Рассмотрим отдельные составляющие комбинированных
моделей.
Застойные зоны. На практике встречаются два вида
застойных зон: застойные зоны без обмена с основным потоком
– «мертвые» зоны − и зоны с обменом между ними и основным
потоком (рис. 3.7).
Рис. 3.7. Схема аппарата с застойными зонами
При наличии обмена между проточной и застойной
зонами возникает задача определения не только объема
застойной зоны, но и эффективности обмена между проточной и
застойной зонами. Характерным признаком наличия в аппарате
застойных зон являются длинные «хвосты» у функции
распределения (на кривой отклика). «Мертвые» застойные зоны
легко определяются индикаторными методами из соотношения
τu =
где
τu
−
среднее
индикаторным методом,
пребывания.
∫ τC(τ)dτ ≠ V
∫ C(τ)dτ V
aп
•
время
τ
=τ,
пребывания,
определенное
− фиктивное среднее время
69
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Фиктивное среднее время пребывания в аппарате можно
представить как
Va
•
V
=
Vпрот
•
+
V
Vзаст
•
= τu +
V
Vзаст
•
V
и
•
•
Vзаст = Vaп − τu V = V(τ − τu ) ,
откуда
τu <
Vaп
•
,
V
•
где V − объемный расход потока; Vaп − объем аппарата, Vпрот
− объем проточной зоны аппарата, Vзаст − объем застойной
зоны аппарата.
Байпасирование (проскок). Пусть требуется определить
долю байпасирующего потока по экспериментальным функциям
отклика.
В случае когда индикатор не попадает в байпасирующий
поток (рис. 3.8), среднее время пребывания индикатора в
аппарате составляет
τu =
∫ τC(τ)dτ = V
∫ C(τ)dτ V
•
или
τu =
Vaп
•
(1 − a) V
,
•
где a =
V1
•
– доля байпасирующего потока.
V
70
aп
2
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Рис. 3.8. Определение байпасирующего
потока: индикатор не попадает в байпасирующий поток
В случае байпасирующего
пребывания определится как
τ=
Vaп
•
потока
среднее
время
,
V
откуда
(3.53)
τ
= 1 − a,
τu
a = 1−
τ
.
τu
Величины τ и τu определяются экспериментально, и с
помощью
соотношения
(3.53)
вычисляется
доля
байпасирующего потока а.
В случае когда индикатор попадает в байпасирующий
поток (рис. 3.9), часть индикатора, прошедшая с
байпасирующим потоком, выйдет из аппарата раньше
индикатора, попавшего в аппарат.
Рис. 3.9. Определение байпасирующего потока:
индикатор
проходит через байпасирующий поток
71
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Если структура потока в аппарате далека от идеального
смешения, то кривая отклика будет иметь два характерных пика
(рис. 3.10). Первый пик соответствует байпасирующему потоку,
второй – реальному аппарату. Определить в этом случае долю
байпасирующего потока можно по следующим соотношениям.
Рис. 3.10. Кривая отклика системы при нанесении
импульсного возмущения для схемы, изображенной на рис. 3.9
Масса индикатора, попавшего в байпасирующий поток
М1:
•
τ1
M1 = V ∫ C(τ)dτ ,
0
масса индикатора M2, попавшего в аппарат, определяется как
• ∞
M 2 = V ∫ C(τ)dτ .
τ1
Принимая,
что
индикатор
разделился
между
байпасирующим потоком и аппаратом пропорционально
соответствующим расходам, получим:
72
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
•
M1 =
M V1
M2 =
M V2
•
= Ma,
V
•
•
= M(1 − a),
V
где М=М1+М2 – общая масса введенного индикатора.
Откуда легко определяется доля байпасирующего
потока.
В случае если в аппарате происходит идеальное
смешение и индикатор, попавший в аппарат мгновенно
появляется на выходе, долю байпасирущего потока можно
определить по решению для соответствующей модели.
Уравнение МИС имеет следующий вид:
dC ' C '− Cвх
+
= 0,
dτ
τ'
здесь С’ – концентрация индикатора на выходе из аппарата; τ ' −
время пребывания индикатора, попавшего в аппарат:
τ' =
•
Vап
=
•
V2
Vап V
•
•
= τ /(1 − a) .
V V2
(3.54)
Решением уравнения (3.54) является:
( τ ')
С ' = Cвх + const exp -τ
Так как индикатор вводится импульсно, Свх = 0, а
const = C '(0) =
M2 M
= (1 − a) = С0 (1 − a) ,
V
V
где C0 − концентрация индикатора в отсутствии байпаса.
Тогда получим
( τ ') = C (1 − a) exp ( −τ(1 − a) τ )
C ' = C0 (1 − a) exp −τ
0
.
(3.55)
Из уравнения материального баланса в точке Z запишем
•
•
V C = C ' V2 ,
73
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
откуда
•
C' =
CV
=
•
V2
C
.
(1 − a)
(3.56)
Подставляя (3.56) в (3.55), получим
 −τ(1 − a) 
С = С0 (1 − a) 2 exp 
.
τ


(3.57)
Сравнение кривых отклика с потоков байпасированием и
без него приведено на рис. 3.11.
По уравнению (3.57) можно определить долю
байпасирующего потока, например, в полулогарифмических
координатах ln
C
= f (τ) , (1−а) будет соответствовать тангенсу
C0
угла наклона этой зависимости.
Рис. 3.11. Кривая отклика системы при нанесении
импульсного возмущения 1 – в отсутствии байпаса;
2 – с байпасированием потока
Комбинированные модели, составленные из параллельно
соединенных зон. Рассмотрим параллельное соединение зон
идеального смешения и идеального вытеснения (рис. 3.12).
Из условия материального баланса в точке Z получаем:
74
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
•
•
•
V1 C1 + V 2 C 2 = V C ,
поэтому концентрация на выходе составляет
•
C = C1
(3.58)
V1
•
V
•
+ C2
V2
•
.
V
Из уравнения (3.58) следует, что такой отклик является
суммой откликов модели идеального смешения и идеального
вытеснения.
Рис. 3.12. Параллельное соединение зон идеального
смешения и идеального вытеснения
На рис. 3.13 изображена кривая отклика на стандартное
возмущение системы, состоящей из параллельного соединения
зон идеального смешения и идеального вытеснения.
Рис. 3.13. Отклик системы, состоящей из параллельно
соединенных зон идеального смешения и идеального
вытеснения на импульсное возмущение
75
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Рассмотренные модели используются для расчета
аппаратов химической технологии в зависимости от типа
аппарата, процесса, протекающего в нем и поставленной цели
исследования. В табл. 3.1 приведены области применения
различных моделей.
Таблица 3.1 Ориентировочные области применения
различных моделей структуры потока в аппарате
1 Модель идеального
вытеснения
2 Модель идеального
смешения
3 Ячеечная модель
4 Ячеечная модель с
обратными
потоками
5 Диффузионная
модель
Трубчатые аппараты с отношением
L / D ≥ 20
Цилиндрические
аппараты
со
сферическим дном в условиях
интенсивного
перемешивания,
аппараты
с
отражательными
перегородками,
барботажные
аппараты близкими соотношениями
L и D
Каскады реакторов с мешалками,
тарельчатые колонны, аппараты с
псевдоожиженным слоем
Тарельчатые и секционированные
насадочные
аппараты,
где
наблюдается заброс вещества в
сторону,
противоположную
направлению
оси
потока
(пульсационные аппараты)
Трубчатые
аппараты,
аппараты
колонного типа с насадкой и без
насадки при осевом рассеивании
вещества
3.2.2. Моделирование теплообменных процессов
Важной задачей при расчете теплообменных аппаратов
является определение поля температуры теплоносителей
Т(τ,x,y,z). Математическая формулировка задач теплообмена
базируется на законах переноса и законах сохранения.
Соответствующие краевые условия определяют начальное
76
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
состояние исследуемого объекта и его взаимодействие с
окружающей средой.
Распределение температуры в потоке жидкости в
различных аппаратах, возникающее вследствие ее движения,
может быть адекватно описано с помощью ранее рассмотренных
моделей движения потоков.
Рассмотрим задачу расчета температур теплоносителей в
аппарате типа «труба в трубе» (схема теплообменного аппарата
приведена на рис. 3.14).
Допустим, что теплоносители движутся по МИВ
прямотоком со скоростью υ , плотность ρ и теплоемкость cр
постоянны. Так как в этом случае температура изменяется
только в продольном направлении, за элементарный объем
можно принять dV = Sdx , где S − площадь поперечного
сечения потока.
Рис. 3.14. Схема теплообменного аппарата
Запишем уравнение теплового баланса для первого
горячего теплоносителя:
υ1T1ρ1c p1S1 − υ1 (T1 + dT1 )ρ1c p1S1 − q12dF = 0 ,
где F – поверхность теплопередачи.
Разделим все члены это уравнения на dV1 = S1dx , а
поток тепла от горячего теплоносителя к холодному q12 выразим
уравнением теплопередачи:
q12 = K(T1 − T2 ) ,
где К – коэффициент теплопередачи.
В итоге получим
υ1
dT1
K(T1 − T2 ) dF
=−
.
ρ1c p1 S1dx
dx
77
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Если поверхность
расстоянием, то
теплообмена
не
меняется
с
dF ПL
=
= П,
dx
L
где П – периметр внутренней трубы; L – длина трубы.
Записывая аналогичное уравнение для второго
теплоносителя, получим следующую систему уравнений:
K(T1 − T2 )
 dT1
П
 dx = − G c
1 p1

,

dT
K(T
−
T
)
2
1
2

=
П
 dx
G 2 c p2
(3.59)
где G = υρS − массовый расход.
Граничные значения для системы уравнений (3.59)
обычно задаются с одной стороны теплообменника: при x=0;
T1 (0) = T1н , T2 (0) = T2н . Тогда решение системы (3.59)
представляет задачу Коши, реализация которой численными
методами является менее трудоемкой, чем задача на граничные
значения, возникающая, например, при противоточном
движении теплоносителей.
Пусть второй теплоноситель движется в обратном
показанному на рис. 3.14, направлении. Уравнение теплового
баланса для этого теплоносителя будет иметь вид
−υ2 T2ρ2 c p2S2 + υ2 (T2 + dT2 )ρ2 c p2S2 + j12dF = 0 .
Таким образом, в случае противоточного движения
теплоносителей математическая модель будет отличаться от
(3.59) знаком перед правой частью второго уравнения:
K(T1 − T2 )
 dT1
П
 dx = − G c
1
p1

.

dT
K(T
T
)
−
2
1
2

П
=−
 dx
G 2c p2
(3.60)
78
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Граничные значения записываются в следующем виде:
при x=0 T1 (0) = T1н ; при x=L T2 (L) = T2н , L – длина
теплообменника.
Алгоритмы
численного
решения
систем
дифференциальных уравнений построены на одновременном
пошаговом вычислении искомых функций вдоль одного
направления по оси абсцисс, поэтому задание граничных
значений на разных концах аппарата при численном решении
системы уравнений (3.60) требует использования метода
итераций. Например, задаваясь последовательно новыми
приближениями T2 (0) = T2к многократно, решается система
дифференциальных уравнений (3.60) до тех пор, пока не будет с
определенной точностью выполнено условие T2 (L) = T2н .
Часто бывает удобно математическое описание
представлять в безразмерных координатах X = x / L :
K(T1 − T2 )
 dT1
=
−
F
 dX
G
c
1
p1

.

−
dT
K(T
T
)
2
1
2

=
F
 dX
G 2 c p2
Рассмотрим
теперь
случай,
когда
движение
теплоносителей описывается диффузионной моделью. Тогда в
уравнении теплового баланса для элементарного объема
необходимо учесть поток тепла, вызванного обратным
перемешиванием,
q L = −DL
dT
ρc P .
dx
Направление
потока
определяется знаком производной. В случае прямотока
dT
υ1T1ρ1cp1S1 − DL1 1 S1ρ1c p1 − υ1 (T1 + dT1 )ρ1c p1S1 +
dx
.
d(T1 + dT1 )
+ D L1
S1ρ1c p1 = K(T1 − T2 )dF.
dx
Преобразуем это уравнение и запишем аналогичное для
второго теплоносителя:
79
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
 dT1
d 2 T1
(T − T )П
υ
=
−K 1 2
D
L1
 1
2
ρ1c p1S1
dx
 dx

2
υ dT2 = D d T2 + K (T1 − T2 )П .
L2
 2 dx
ρ2 c p2S2
dx 2

(3.61)
Или для безразмерной координаты:
 dT1
(T − T )F
1 d 2 T1
υ
=
−K 1 2
 1
2
cp1G1
 dX Pe1 dX

2
υ dT2 = 1 d T2 + K (T1 − T2 )F .
 2 dX Pe2 dX 2
c p2 G 2

Так как уравнения диффузионной модели являются
дифференциальными уравнениями второго порядка, для
решения системы (3.61) необходимо записать по 2 граничных
условия для каждого уравнения. Обычно эти граничные условия
имеют следующий вид (по Данквертсу):
для теплоносителя 1:
при х=0 υ1 (T1H − T1 ) = − D L1
dT1
dx
; при x=L
x =0
dT1
dx
=0;
x =L
для теплоносителя 2:
при х=0 υ2 (T2H − T2 ) = − D L2
dT2
dx
; при x=L
x =0
dT2
dx
=0.
x =L
Пояснения к граничным условиям данного типа
приводятся на рис. 3.15. Таким образом, решение уравнения
диффузионной модели связано с решением задачи на граничные
значения и выполняется итерационно численными методами.
Особенностью получаемого решения с граничными значениями
по Данквертcу является скачкообразное изменение температуры
теплоносителя при его входе в аппарат, например для первого
теплоносителя − от T1H до T1 .
80
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Рис. 3.15. К постановке граничных условий
В случае противоточного движения теплоносителей
согласно диффузионной модели система уравнений будет иметь
вид
 dT1
d 2 T1
(T − T )П
υ
=
D
−K 1 2
L1
 1
2
dx
ρ1c p1S1
 dx
.

2
dT
d
T
(T
−
T
)
П
υ
2
2
= − D L2
−K 1 2
2
 2 dx
ρ2 c p2S2
dx

А граничные условия:
- для теплоносителя 1:
при х=0 υ1 (T1H − T1 ) = − D L1
dT1
dx
dT1
dx
x =0
dT2
dx
x =L
; при x=L
=0;
x =L
- для теплоносителя 2:
при х=L υ2 (T2H − T2 ) = D L2
dT2
dx
; при x=0
=0;
x =0
Структура потоков, близкая к модели идеального
смешения, реализуется в аппаратах с интенсивным
перемешиванием, например в аппарате с мешалкой,
изображенном на рис. 3.16. Предполагая однородное поле
81
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
температур по всему объему аппарата, получим следующее
уравнение МИС:
(T1 − T1H ) =
K(T1 − TП )F
,
G1cp1
где ТП – температура пара, подаваемого в рубашку аппарата.
Рис. 3.16. Схема аппарата с мешалкой
Далее несложно получить описание теплообмена по
ячеечной модели. Математическое описание включает
уравнения теплового баланса для каждой ячейки, например для
i-ой ячейки
Gc p (T1, j − T1, j−1 ) = K j
F
(TП − T1, j ),
m
j=1,m.
Структура потока оказывает существенное влияние на
распределение температуры теплоносителей по длине
теплообменника и соответственно на движущую силу процесса.
На рис. 3.17 представлены профили температуры в ядре потока
вдоль поверхности теплоотдачи для различных моделей
структуры потока в аппарате одного из теплоносителей.
Температура второго теплоносителя была неизменной и
равнялась ТП, что имеет место в случае теплообмена при
конденсации паров. Видно, что наибольшая движущая сила и
наибольшее
изменение
температуры
теплоносителя
соответствуют МИВ, а наименьшее − МИС. Диффузионная и
ячеечная модель занимают промежуточные положения. Как
будет показано далее, подобные выводы верны для многих
процессов массообмена и химического превращения.
82
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Рис. 3.17. Профили температуры в ядре потока
вдоль поверхности теплоотдачи для
различных моделей структуры потока в аппарате
Полученных ранее моделей оказывается недостаточно
для решения задач теплообмена, так как входящие в них
теплофизические свойства теплоносителей и коэффициенты
переноса тепла оказываются неопределенными. Причем их
значения будут сильно влиять на получаемые результаты
расчета. Для определения этих характеристик можно
использовать как экспериментальные данные, так и различные
модели. Последнее удобнее, так как позволяет учитывать
изменение теплофизических и переносных свойств с
изменением условий протекания процесса.
Напомним, что коэффициент теплопередачи К
выражается через коэффициенты теплоотдачи α1, α2,
термического сопротивления стенки и загрязнений.
Для плоской стенки коэффициент теплопередачи
определяется как
K=
1
δ 1
1
+∑ +
α1
λ α2
,
где α1 – коэффициент теплоотдачи от горячего теплоносителя; δ
– толщина теплопередающей стенки аппарата, м; λ −
коэффициент теплопроводности материала стенки; α2 −
83
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
коэффициент
теплоотдачи
от
стенки
к
холодному
теплоносителю.
Коэффициент теплоотдачи является сложной функцией
от многих факторов:
• скорости жидкости, ее плотности и вязкости, то есть
переменных, определяющих режим течения жидкости;
• теплофизических свойств жидкости (удельной
теплоемкости, теплопроводности, а также коэффициента
объемного расширения;
• геометрических параметров – формы и определяющих
размеров стенки (для труб – их диаметр и длина), а также
шероховатости стенки.
Практически коэффициенты теплоотдачи определяют по
критериальным уравнениям, связывающим критерий Нуссельта
Nu =
αl
, с другими критериями подобия, значения которых
λ
влияют на процесс. Вид этих уравнений зависит от типа
процесса теплоотдачи (при обтекании пучка труб, при
конденсации, при кипении и т.д.).
Многие из факторов, влияющих на коэффициент
теплоотдачи, сами зависят от температуры, поэтому часто в
расчетах теплообменных процессов приходится учитывать
зависимость коэффициента теплоотдачи от температуры.
3.2.3. Моделирование массообменных процессов
В качестве аппаратов для массообменных процессов
(абсорбция, ректификация, экстракция и др.) широко
используются тарельчатые и насадочные колонны. Характер
математического моделирования массообменных процессов в
аппаратах колонного типа весьма сходен. В связи с этим
основные принципы моделирования будут продемонстрированы
на примерах насадочной абсорбционной колонны и тарельчатой
колонны для ректификации двухкомпонентной (бинарной)
смеси.
Особенностью протекания массообменных процессов
является возможное изменение расходов фаз по высоте колонны
в результате перераспределения компонентов между ними. При
84
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
составлении материального баланса для элементарного объема
удобно за расходы фаз принять величины, не изменяющиеся по
высоте колонны. В противном случае расходы либо попадут под
оператор дифференцирования, либо необходимо будет
записывать дополнительные уравнения, что усложнит
математическую модель. Выбранные за расходы фаз величины в
дальнейшем определяют способ выражения состава смесей и
запись законов переноса и равновесия. Например, если за расход
фазы принят расход всей смеси, выраженный в кмоль/с или кг/с,
то состав фазы необходимо определять в мольных или массовых
долях соответственно. При этом нужно соблюдать соответствие
в размерностях коэффициентов переноса (массоотдачи,
массопередачи, диффузии) и движущей силы процесса.
Рассмотрим процесс физической абсорбции при условии,
что потоки газа и жидкости движутся согласно МИВ. Выделим
элементарный объем dV = Sdz. Для процесса абсорбции расходы
следует выражать через расходы инертных компонентов, так как
они не изменяются в процессе. Пусть G и L – массовые расходы
инертного компонента в газовой и жидкой фазах
соответственно.
Тогда
концентрации
распределяемого
компонента А в газовой Y и жидкой X фазах необходимо
выражать в относительных массовых концентрациях кгА/кгин.
Произведение GY или LX определяют массовый расход
распределяемого компонента в газовой и жидкой фазах
соответственно.
Для МИВ запишем уравнение материального баланса по
распределяемому компоненту для элемента объемом Sdz (рис.
3.18), где S – площадь поперечного сечения колонны:
GY − G(Y + dY) − jгyx dF = 0 − газовая фаза,
−LX + L(X + dX) + jгyx dF = 0 − жидкая фаза,
jгyx − межфазный поток массы, который определим уравнением
массопередачи
jгyx = K Y (Y − Y* (X)) ;
KY
− коэффициент
массопередачи; dF − поверхность массопередачи. Поделив
уравнения материального баланса на величину элементарного
объема, получим
85
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
 G dY
*
 S dz = − K VY (Y − Y (X))

,
 L dX = − K (Y − Y* (X))
VY
 S dz
(3.62)
где K VY = K Y a
a=
− объемный коэффициент массопередачи;
dF
− удельная поверхность контакта фаз, которая зависит
dV
от
типа
контактного
экспериментально.
устройства
и
определяется
Рис. 3.18. Структура потоков фаз в абсорбере
Граничные условия для системы уравнений (3.62)
записываются в виде
Y(0) = Yн
X(H) = X н
.
При расчете коэффициента массопередачи необходимо
соблюсти соответствие его размерности с размерностью
движущей силы в математической модели (3.62). Коэффициент
массопередачи определяется согласно соотношению
KY =
1
m 1
+
βX βY
,
где β − коэффициент массоотдачи, который, как и коэффициент
теплоотдачи, является сложной функцией режима движения и
физико-химических свойств фаз, геометрии аппарата; m –
коэффициент распределения. Коэффициенты массоотдачи
86
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
обычно рассчитываются по критериальным уравнениям, причем
часто
размерность
рассчитываемых
коэффициентов
массоотдачи (м/с), соответствует движущей силе, выраженной в
объемно массовых концентрациях (кг/м3), поэтому перед
подстановкой в (3.62) определенные таким образом
коэффициенты необходимо пересчитать в соответствующую
размерность.
Описание процесса абсорбции (как и любого
массообменного процесса) включает в себя описание фазового
равновесия y* = f(x). Без уравнения равновесия, которое
определяет связь концентраций распределяемого компонента в
фазах, невозможно решить систему уравнений (3.62)Точность
описания условий фазового равновесия в значительной степени
определяет адекватность получаемых по данной модели
результатов. Более того, существуют модели, например модель
теоретической тарелки, результаты расчета по которым целиком
определяются точностью описания фазового равновесия.
Для случая плохой растворимости газа в жидкости
применим закон Генри
y* =
E
x,
P
(3.63)
где Е – константа Генри значение которой определяют по
справочным данным [9,10]; P − общее давление газовой смеси.
Коэффициент распределения в данном случае определяется как
m=
E
. В соотношении (3.63) концентрации выражаются в
P
мольных долях, поэтому перед использованием закона Генри в
математической модели (3.62) его необходимо переписать в
соответствующих концентрациях. Для относительных массовых
концентраций получим
Y* =
Mж Е
Mг Р
X
M
 Е
1 + Х ж + X 1 − 
Mа
 Р
,
где Мж, Мг, Ма – молекулярные массы абсорбента, инертной
газовой фазы и поглощаемого компонента. В случае малых
87
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
концентраций поглощаемого компонента в абсорбенте, это
соотношение можно упростить:
Mж Е
X,
Mг Р
M Е
тогда коэффициент распределения m = ж .
Mг Р
Y* =
Далее рассмотрим описание структуры потоков газовой
и жидкой фаз в абсорбере диффузионной моделью. Поток массы
за счет обратного перемешивания имеет вид
j = −DLX
dC X
.
dz
В
этом
соотношении
концентрация
СX
объемномассовая
размерностью
кгA/м3
и
связана
относительной массовой концентрацией соотношением:
CX = XС ИН ,
где
СИН
−
объемномассовая
концентрация
–
с
абсорбента,
кг ин / м3 .
В случае малых концентраций распределяемого
компонента А, что обычно характерно для растворения газов
жидкостях, СИН ≈ ρX , где ρX − плотность чистого абсорбента
при температуре и давлении в абсорбере. Тогда ρX X ≈ CX и
j = − D LX ρX
dX
.
dz
Для
газовой
фазы
равенство
обемномассовой
концентрации инертного компонента и его плотности в чистом
виде не так очевидно. Однако, учитывая, что в процессе
абсорбции концентрация инертного носителя увеличивается, в
среднем можно принять СИН ≈ ρY , где ρY − плотность чистого
инертного носителя газовой фазы при температуре и давлении в
абсорбере, и записать
j = − D LY ρY
dY
.
dz
Запишем уравнения материального
выделенного элементарного объема (рис. 3.18):
88
баланса
для
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
dY

GY − D LY dz ρYS − G(Y + dY) +

+ D d(Y + dY) ρ S − K (Y − Y *(X))dV = 0
Y
VY
 LY
dz
.

− LX − D dX ρ S + L(X + dX) +
LX
X

dz

d(X + dX)
+ DLX
ρXS + K VY (Y − Y * (X))dV = 0
dz

Проведя несложные
уравнений, получим
преобразования
этой
системы
 G dY
d2Y
=
ρY − K VY (Y − Y * (X))
D
LY
 S dz
dz 2
.

2
 L dX = − D d X ρ − K (Y − Y *(X))
LX
X
VY
 S dz
dz 2
Граничные условия:
− для фазы Y:
при z=0 G(YH − Y) = −SρY D LY
dY
; при z=H
dz z =0
dY
=0;
dz z = H
− для фазы X:
при z=H L(X H − X) = Sρ x DLX
dX
; при z=0
dz z = H
dX
= 0.
dz z =0
Несложно
получить
математические
массообмена для случая идеального смешения в фазах:
G ( YH − YК ) = K VY (Y − Y *(X К ))VАП
,

L ( X К − X Н ) = K VY (Y − Y *(X К ))VАП
89
модели
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
где VАП − объем аппарата.
Согласно ячеечной модели, абсорбер представляют в
виде m последовательно соединенных ячеек идеального
смешения. Для i-ой ячейки:
VАП

G ( Yi −1 − Yi ) = K VYi (Yi − Y *(X i )) m
.

L ( X − X ) = K (Y − Y *(X )) VАП
i
i −1
VYi
i
i

m
Моделирование процесса хемосорбции
При
хемосорбции
(абсорбция,
сопровождаемая
химической реакцией) абсорбируемый компонент связывается в
жидкой фазе в виде химического соединения. Протекание
химической реакции в процессе абсорбции оказывает влияние
на движущую силу процесса массопередачи, поскольку
изменяется концентрация распределяемого вещества в жидкой
фазе.
Как и в случае физической абсорбции, рассматриваемая
система будет состоять из двух фаз − газ и жидкость. В газовой
фазе должно быть минимум два компонента, в жидкой минимум
три (если сам абсорбент вступает в обратимую химическую
реакцию с поглощаемым компонентом). Например, известный
процесс поглощения СО2 из воздуха водным раствором
моноэтаноламина, в котором газовая фаза состоит из двух
компонентов – воздуха и углекислого газа, а жидкая фаза – из
четырех компонентов: моноэтаноламина, воды, СО2 и продукта
реакции.
Рассмотрим абсорбер, в котором протекает реакция
между поглощаемым компонентом А и хемосорбентом В,
растворенным в жидкой фазе:
А + В 
→C .
Сам абсорбент является инертным, то есть в химическую
реакцию не вступает. Кинетика этой химической реакции
определена соотношением r = kCnAA CnBB .
Рассмотрим только случай движения фаз по МИВ. Для
других видов структур потоков переход от физической
90
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
абсорбции к хемосорбции будет аналогичным. Запишем
уравнение материального баланса для жидкой фазы. Для
газовой фазы уравнение материального баланса останется таким
же, как и для физической абсорбции. В жидкой фазе происходят
изменения концентраций компонентов А и В, поэтому
необходимо записать два уравнения материального баланса:
− LX A + L(X A + dX A ) + K VY (YA − Y *(X A ))dV − rdV = 0
,
− LX B + L(X B + dX B ) − rdV = 0 .
Уравнения для нахождения скорости химической
реакции необходимо переписать в относительных массовых
концентрациях. Учитывая переход, показанный ранее
A +nB )
r = kCnAA CnBB = kρ(n
X nAA X nBB = kX nAA X nBB .
x
Здесь необходимо заметить, что в отличие от k, которая
зависит только от температуры, модифицированная константа
скорости k будет зависеть еще и от давления. В результате
получим следующую математическую модель:
L dX A
= − K VY (YA − Y * (X A )) + kX An A X Bn B
S dz
L dX B
= kX nAA X nBB
S dz
G dYA
= − K VY (YA − Y * (X A ))
S dz
с соответствующими граничными значениями:
X A (H) = X Aн
X B (H) = X Bн .
YA (0) = YAн
Если растворимость компонента А мала, то есть XА→0,
тогда можно пренебречь изменением концентрации компонента
В (XВ=const) и не рассматривать второе уравнение.
Неизотермическая абсорбция
При растворении веществ друг в друге часто выделяется
или поглощается определенное количество теплоты, называемой
91
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
теплотой растворения. Если тепловой эффект достигает
достаточно больших значений, то это приводит к значительному
изменению температуры в абсорбере, что сказывается на
растворимости поглощаемого компонента.
Составим математическое описание неизотермической
абсорбции для МИВ. Уравнения для изменения концентраций
компонентов в фазах запишутся как в (3.62). Для учета
изменения температуры фаз составим уравнения теплового
баланса:
− для жидкой фазы:
− LTX c P + L(TX + dTX )c P + K VY (Y − Y * (X, TX ))qdV − Ka(TX − TY )dV = 0
;
− для газовой фазы:
GTY c PY − G(TY + dTY )c PY + Ka(TX − TY )dV = 0 ;
где q − тепловой эффект растворения (может быть разного
знака). Некоторым приближением в этих уравнениях является
то, что, за массовые расходы фаз приняты расходы инертных
носителей. Поделив на элементарный объем, окончательно
получим:
 L dX
 S dz = − K VY (Y − Y *(X, TX ))

 G dY = − K (Y − Y *(X, T ))
VY
X
 S dz
 Lc
 PX dTX = − K (Y − Y *(X, T ))q + Ka(T − T )
VY
X
X
Y
 S dz
 Gc
 PY dTY = Ka(TX − TY )
 S dz
и граничные условия:
X(H) = X н
Y(0) = Yн
TX (H) = TXн
TY (0) = TYн .
92
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
В условии равновесия надо учесть зависимость коэффициента
Генри и давления от температуры: Y* =
M ж Е (TX )
X.
M г Р(TY )
Математическое описание процесса ректификации в
тарельчатой колонне
Рассмотрим тарельчатую ректификационную колонну, в
которой осуществляется ступенчатый контакт фаз, содержащую
N тарелок. На рис. 3.19 представлена соответствующая
ректификационная установка, включающая в себя саму
тарельчатую ректификационную колонну, куб–испаритель
(кипятильник) и дефлегматор.
Кипятильник помещают внизу колонны для получения
парового потока. Жидкость, стекающая с нижней тарелки,
поступает в кипятильник, где испаряется. Образовавшийся пар
поднимается вверх по колонне. Дефлегматор помещают вверху
колонны для конденсации выходящих из нее паров. Он может
быть двух типов: полный и парциальный. При применении
полного дефлегматора весь пар, поднимающийся с верхней
тарелки колонны конденсируется. Часть образовавшейся
жидкости отводится в виде верхнего продукта (дистиллята),
остальная жидкость, называемая флегмой, возвращается в
колонну.
93
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
D, xd
Дефлегматор
Вода
y1
0
x0
1
2
3
.
.
.
.
.
F, xf
N-2
N-1
N
xN
yN+1
Пар
N+1
Конденсат
Куб
W, xw
Рис. 3.19. Схема ректификационной установки
В парциальном дефлегматоре пар, отходящий с верхней
тарелки
колонны,
конденсируется
только
частично,
образовавшийся конденсат делится на дистиллят и флегму,
которая поступает обратно в колонну на орошение.
Несконденсировавшиеся пары направляются на дальнейшую
переработку или утилизацию.
Применительно к тарельчатой ректификационной
колонне наиболее пригоден метод расчета «от тарелки к
тарелке». Согласно этому методу все уравнения процесса
ректификации решаются для каждой ступени разделения
(тарелки) в отдельности. Нумерацию тарелок колонны
традиционно ведут сверху вниз, так как при расчете процесса
ректификации
необходимо
учитывать
дефлегматор
и
кипятильник; дефлегматору соответственно присваивают
нулевой номер, а кипятильнику – N + 1.
При моделировании тарельчатых ректификационных
колонн методом от тарелки к тарелке используются два
варианта алгоритмов:
• однонаправленный расчет (снизу вверх или сверху
вниз);
94
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
• расчет от концов колонны к тарелке питания.
Для описания многокомпонентной ректификации
наиболее предпочтительным является второй вариант. Решение
ищется итерационным способом. Согласно этому алгоритму в
начале каждой итерации заданными являются составы
продуктов разделения xd и xw. С обоих концов колонны для
каждой ступени разделения решается система уравнений,
описывающих процесс массообмена на тарелке. Согласование
решения проводится по составу пара над тарелкой питания,
который получается в расчете сверху и снизу колонны. Если
составы паровой фазы не совпадают, то производится
корректировка составов продуктов разделения xd и xw, и
начинается новый итерационный цикл.
При использовании итерационного метода, как правило,
возникают
трудности
с
обеспечением
сходимости
итерационного процесса, поэтому для обеспечения сходимости
используются специальные методы, смысл которых заключается
в обработке в определенном порядке вычисленных значений
переменных
с
целью
получения
улучшенных
(корректированных) значений для следующего приближения,
которые ближе к истинным величинам. Наиболее широко
распространен θ-метод сходимости.
В данном методе перед новой итерацией производится
корректировка начальных приближений расхода каждого
компонента выходящего из верхней di (мольный расход
компонента i в дистилляте, моль/с]) и нижней wi частей
колонны:
(d i )1 =
Fx f i
,
1 + θ(w i / d i )0
w 
( wi )1 = θ  i  ( d i )1 .
 d i 0
Поправочный множитель θ определяется из решения
уравнения
m
Fx f i
−D,
i =1 1 + θ(w i / d i ) 0
g(θ) = ∑
95
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
где F и D – расходы питания и дистиллята соответственно,
моль/с; xfi – мольная концентрация компонента i в питании.
Искомым значением θ является положительный корень, при
котором g (θ) = 0 .
Скорректированные значения мольных расходов di и wi
переводятся в мольные концентрации xdi и xwi. Концентрации
потоков пара и жидкости, поступающих в колонну (yN+1, x0) и
выходящих из нее (y1, xN), зависят от типа выбранного
дефлегматора и кипятильника. Для полного дефлегматора
состав отходящего с верхней тарелки пара равен составу
дистиллята и поступающей в колонну на орошение флегмы
(y1=x0=xd). Во втором случае пар в дефлегматоре y0 находится в
равновесии с жидкостью, поэтому его рассматривают как
дополнительную равновесную ступень (y0 = y*(x0), y0 = xd). В
производственной практике больше распространен первый
вариант.
В случае парциального кипятильника пар, отходящий из
него, находится в равновесии с уходящей из него жидкостью,
соответственно его также можно рассчитывать как равновесную
ступень (N+1 тарелку). В этом случае концентрация пара
определяется как равновесная с уходящей кубовой жидкостью
yN+1=y*(xw), а концентрация xN определяется соответственно из
материального баланса. Однако в расчетах часто используется
кипятильник полного испарения, в котором состав выходящего
парового потока равен составу поступающей в него жидкости.
Следовательно, в этом случае yN+1=xN=xw. В расчетах для
упрощения рассматривают кипятильник полного испарения.
Для описания процесса массообмена на тарелке по
жидкой фазе часто используется ячеечная модель структуры
потока (рис. 3.20), при этом тарелка делится по ходу жидкости
на
ряд
ячеек
m.
Как
показали
многочисленные
экспериментальные исследования, для описания движения
парового потока в массообменном пространстве тарелки
адекватной является модель идеального вытеснения.
96
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Рис. 3.20. Схема потоков пара и жидкости
на однопоточной тарелке.
Ячеечная структура потока жидкости на тарелке:
L – поток жидкости; G – поток пара;
α – номер тарелки; j – номер ячейки
Тогда для определения изменений концентраций в j
ячейке тарелки α запишем следующую систему уравнений:
L(x j+1 − x j ) = G(y αj (h) − y αj (0))
G dyαj
= − K y a ( yα , j − y * (x αj ) ) .
mS dz
При этом предполагается, что:
• расходы жидкости и пара для каждой ячейки
постоянны;
• на входе в каждую ячейку пар имеет один и тот же
y α+1, j (0) = y α+1
.
На выходе из ячеек состав пара по ходу движения
жидкости на тарелке будет меняться, поэтому после каждой
ступени разделения производится усреднение состава пара:
состав
m
∑y
yα =
j=1
m
αj
.
Для описания ступени разделения необходимо
определить число ячеек полного перемешивания, на которое
условно
можно
разделить
тарелку.
Число
ячеек
преимущественно зависит от гидродинамической обстановки.
97
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Скорость газовой фазы в интервале 0.4 − 0.97 м/с не оказывает
заметного влияния на интенсивность перемешивания жидкой
фазы. Наибольшее влияние на интенсивность перемешивания
жидкости оказывают расход жидкости и высота сливной
перегородки; с увеличением расхода жидкости перемешивание
ослабевает, с увеличением высоты сливной перегородки −
интенсифицируется. Число ячеек полного перемешивания,
например, для ситчатой тарелки рассчитывается по уравнению
[11]:
m = 0, 45 + lT (8, 4 − 0.036h с.п. ) − 0, 0028G / L,
где lT – длина пути жидкости на тарелке, м; hс.п. – высота
сливной перегородки, м. Существуют и другие выражения для
определения числа ячеек полного перемешивания n в колоннах с
ситчатыми тарелками диаметром до 600 мм [12].
Для колонн диаметром более 600 мм с ситчатыми,
колпачковыми и клапанными тарелками отсутствуют надежные
данные по продольному перемешиванию жидкости, поэтому с
достаточной степенью приближения можно считать, что одна
ячейка перемешивания соответствует длине пути жидкости l =
300 − 400 мм. Обычно принимают l = 350 мм, тогда число ячеек
полного перемешивания определяется как отношение длины
пути жидкости на тарелке lТ к длине l. Длина пути жидкости lТ
определяется как расстояние между переливными устройствами.
В случае установки в колонне многопоточных тарелок
длина пути жидкости lТ сокращается, что приводит
соответственно к сокращению числа ячеек. И при некотором
диаметре колонны длина пути жидкости на тарелке lТ может
быть равной длине пути жидкости l для ячейки перемешивания.
В этом случае для описания ступени разделения достаточно
одной ячейки идеального перемешивания, что существенно
сокращает размерность задачи.
Выше
были
рассмотрены
различные
примеры
составления математических моделей. Можно построить
общую, исходную схему (рис. 3.21).
98
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Рис. 3.21. Общая схема составления математических
моделей
3.2.4. Математические модели химических реакторов
Химические превращения являются ядром любой
химической технологии, поэтому эффективность проведение
этих процессов существенно влияет на качество конечного
продукта и экономические показатели всей технологической
цепочки. Изучением закономерностей протекания химических
реакций занимаются специалисты различных областей химии.
Однако данных, получаемых этими специалистами в виде
механизмов реакций, уравнений формальной химической
кинетики и численных значений параметров в этих уравнениях,
оказывается недостаточно для проектирования эффективных
промышленных способов проведения химических процессов.
Для успешной реализации научной идеи, полученной и
исследованной в химической лаборатории, в промышленном
масштабе необходимо уметь учитывать влияние на протекание
химических реакций таких факторов, как гидродинамическая
структура потока, тепло- и массообмен. Степень влияния этих
факторов может быть разной в зависимости от типа реакций и
условий их протекания. Например, массообмен оказывает
существенное влияние в основном на протекание гетерогенных
реакций и может не учитываться при рассмотрении реакций,
протекающих гомогенных средах. Здесь будут рассмотрены
некоторые закономерности влияния структуры потоков тепло- и
массообмена на протекание простых и сложных реакций.
В зависимости от тепловых условий протекания
химических реакций различают следующие типы реакторов:
99
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
− изотермический – реакция происходит при постоянной
температуре;
− адиабатический – отсутствует теплообмен с
окружающей средой;
− политропический – наличие теплообмена с
окружающей средой.
Рассмотрим стационарные модели химических реакторов
с различной структурой потока на примере изотермического
реактора, когда при помощи подвода или отвода тепла в
реакторе поддерживают постоянную температуру в течение
всего процесса. Пусть в реакторе протекает простая гомогенная
реакция A → B .
В случае идеального вытеснения материальный баланс
для элементарного объема реактора имеет вид (рис.3.22)
SυCA − Sυ(C A + dCA ) − rdV = 0 .
Рис.3.22. К построению моделей химических реакторов
Принимая, что dV = Sdx и r = kC A , после несложных
преобразований получим
υ
dC A
= − kCA
dx
(3.64)
и граничное условие CA (0) = CA Н .
Часто уравнение (3.64) записывают не в неподвижной
системе координат, связанной с аппаратом, а в системе, которая
движется с постоянной скоростью υ . Тогда уравнение (3.64)
преобразуется к виду
100
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
dC A
= −kC A .
dτ
(3.65)
Между переменными уравнений
существует очевидная связь: υτ = x .
Для диффузионной модели:
(3.64)
и
(3.65)
dCA
d(C A + dCA )
− Sυ(CA + dC A ) + SDL
− rdV = 0;
dx
dx
dC
d 2 CA
υ A = DL
− kCA .
dx
dx 2
SυCA − SD L
Граничные условия задаются аналогично (3.39).
Для модели идеального смешения получим
C A − CA Н
τ
= − kC A .
(3.66)
Теперь рассмотрим другие типы реакторов только на
примере модели идеального вытеснения. Для других вариантов
структур потоков эти обобщения можно провести аналогичным
образом.
При адиабатическом режиме в реакторе отсутствует
теплообмен с окружающей средой, и тепло химической реакции
полностью расходуется на изменение температуры реакционной
смеси. Запишем тепловой баланс для элементарного объема
реактора:
Sυρc p T − Sυρc p (T + dT) + rq r dV = 0 ,
где q r − удельная теплота химической реакции (может иметь
различный знак).
Учитывая, что уравнение материального баланса будет
аналогично (3.65), получим следующую систему уравнений:
dC A
= − k(T)CA
dτ
dT
υρcp
= k(T)C A q r
dx
(3.67)
101
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
с граничными условиями: CA (0) = CA H и T(0) = TН . Систему
уравнений (3.67) необходимо дополнить уравнением Аррениуса
(3.33), устанавливающим связь между температурой и
константой скорости химической реакции.
При политропическом режиме температура в реакторе
также непостоянна, но при этом часть тепла может отводиться
от реакционной смеси или подводиться к ней. Тогда в
уравнении теплового баланса для элементарного объема надо
учесть поток тепла теплопередачей. В результате получим:
dC A
= − k(T)CA
dτ
dT
υρcp
= k(T)CA q r − K(T − T0 )П,
dx
где T0 − температура среды с внешней стороны стенки
реактора; K − коэффициент теплопередачи; П − поверхность
теплопередачи.
3.2.5. Влияние структуры потоков на протекание
химических реакций
Важными критериями эффективности проведения
химических превращений в промышленных аппаратах являются
степень превращения (конверсия) и селективность.
Степень превращения Х какого-либо реагента А равна
доле превратившегося в продукты вещества от общего
начального количества этого вещества. Например, для реакции
A+B→C
Степень превращения вещества А можно выразить
следующим образом:
XA =
M 0A − M A
M
= 1 − 0A ,
0
MA
MA
а степень превращения вещества В так:
XB = 1 −
102
MB
,
M 0B
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
где МА, МВ – количество вещества А и В в текущий момент,
M 0A , M B0 − количество вещества А и В в начальный момент.
Селективность Sе − выход целевого продукта на
затраченное исходное вещество:
Se =
Количество реагента, перешедшее в целевой продукт
.
Общее количество превратившегося реагента
Таким образом, величина Se показывает, какая доля
превратившегося исходного вещества затрачена на основную
реакцию – на образование целевого продукта.
Вначале рассмотрим влияние структуры потоков в
реакторе на протекание простых реакций типа A → B .
Анализировать в этом случае имеет смысл только степень
превращения.
Для модели идеального вытеснения (3.65) решение будет
иметь вид
x
CA (τ) = C0A exp(−kτ) или CA (x) = C0A exp(−k ) ,
υ
(3.68)
где C0A − концентрация компонента А на входе в аппарат.
Степень превращения определится как
X A (τ) = 1 − exp(− kτ) .
Для идеального смешения (3.66) соответственно получим
CA ( τ ) =
C0A
1+ kτ
и X( τ ) =
kτ
1+ kτ
(3.69)
На рис. 3.23 представлены поведение концентраций
реагента в реакторе идеального смешения и реакторе
идеального вытеснения. Если рассмотреть одно и то же
произвольное поперечное сечение обоих аппаратов, то значение
CA в аппарате идеального вытеснения окажется больше
соответствующего значения в аппарате смешения, поэтому
скорость реакции, пропорциональная CA , во всех сечениях
аппарата, кроме последнего (на выходе), также будет больше в
103
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
аппарате
идеального
вытеснения.
Получить
картину,
показанную на рис. 3.23, где начальные и конечные
концентрации одинаковы для обоих аппаратов, можно только в
том случае, если объем аппарата смешения больше объема
аппарата вытеснения.
CA
CA0
1
CAK
2
L
x
Рис. 3.23. Изменение концентрации реагента
по длине аппарата:
1 – идеального вытеснения; 2 – идеального смешения
Количественные оценки эффективности разных структур
потоков можно провести, задав ряд значений kτ и рассчитав
степень превращения в потоках идеального вытеснения Хвыт и
идеального смешения Хсм. Результаты такого расчета
представлены в табл. 3.2.
Таблица 3.2. Степень превращения в реакторах
идеального смешения и идеального вытеснения
kτ
Xвыт
Xсм
0,5
0,394
0,333
1
0,632
0,5
2
0,865
0,667
104
4
0,982
0,800
6
0,9975
0,856
10
0,99996
0,909
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
При малых значениях kτ , соответствующих случаям
малого объема аппарата или медленной реакции, разница не
очень велика, однако с ростом kτ она может оказаться очень
большой. Картина становится наглядной при сравнении
результатов проектного расчета, когда задается величина
степени превращения и по времени пребывания рассчитываются
необходимые объемы реакторов (табл. 3.3).
Таблица 3.3. Объемы реакторов идеального смешения и
идеального вытеснения в зависимости от степени превращения
X
Vсм/ Vвыт
0,5
1,5
0,7
2
0,9
4
0,95
6
0,99
22
0,999
140
Видно, что при конверсиях близких к единице, объемы
аппаратов идеального вытеснения и смешения могут отличаться
друг от друга более чем в сто раз. Таким образом, важно понять,
что на протекание реакций в промышленных аппаратах
существенное влияние оказывает структура потока, при этом
аппарат идеального вытеснения обеспечивает большую
эффективность процесса. Вывод, что в аппарате идеального
вытеснения глубина превращения выше, чем в аппарате
идеального смешения, справедлив для изотермических
необратимых и обратимых реакций любого порядка (кроме
нулевого), а также для большинства тепло- и массообменных
процессов. Последнее утверждение продемонстрировано на рис.
реальных
потоков,
описываемых
3.17.
Структуры
диффузионной или ячеечной моделью, будут занимать
промежуточное по эффективности положение между идеальным
вытеснением и смешением.
Теперь рассмотрим сложные многостадийные реакции.
Характер влияния структуры потока на сложные реакции
отличается большим многообразием. При проведении сложных
реакций с побочными стадиями одной из главных задач
является достижение высокой селективности. Зачастую ради
этого
жертвуют
степенью
превращения,
так
как
непрореагировавшие компоненты можно выделить из
вышедшей из аппарата смеси и вернуть в начало процесса.
Тогда дополнительные затраты будут связаны только с
дополнительными процессами разделения. Низкая же
105
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
селективность будет означать, что часть ценных реагентов
расходуется впустую, необходимо выделять побочные продукты
из получаемой в реакторе смеси, использовать процессы
утилизации побочных продуктов.
Рассмотрим два случая реакций с последовательной и
параллельной побочной стадией. Пусть последовательная
реакция имеет следующий вид:
1 по А; k1
A 
→B
,

1 по B; k 2
→C
B 
где В − целевой продукт; С – отброс. Над стрелками обозначен
порядок каждой реакции и соответствующие константы
скорости. Примем условие k1 = k 2 = k . Тогда для идеального
вытеснения система уравнений кинетики химических реакций
будет иметь вид
 dC A
 dτ = − kC A
.

dC
B

= kC A − kCB
 dτ
(3.70)
Концентрация компонента С определяется из уравнения
материального баланса CC (τ) = CA0 − C A (τ) − C B (τ) .
Тогда получаем следующее решение:
CA (τ) = C0A exp(−kτ);
CB (τ) = CA0 kτ exp(−kτ);
CC (τ) = CA0 [1 − (1 + kτ) exp(−kτ) ] ,
(3.71)
Откуда для степени превращения компонента А и
селективности легко получить уравнения
X A (τ) = 1 − exp(− kτ);
S(τ) =
CB (τ)
kτ
=
.
CB (τ) + CC (τ) 1 − exp(− kτ)
(3.72)
106
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Для идеального смешения уравнения для изменения
концентраций компонента А и В будут иметь следующий вид:
C A − C0A
= − kCA
τ
C B − C0B
= kCA − kCB .
τ
(3.73)
Для компонента С верно записанное ранее уравнение
материального баланса. Легко получить решение этих
уравнений:
C0A
CA ( τ ) =
;
1 + kτ
C0A k τ
CB ( τ ) =
;
(1 + k τ ) 2
CC ( τ ) =
C0A (k τ ) 2
.
(1 + k τ ) 2
(3.74)
Откуда
kτ
1+ kτ
1
Se( τ ) =
= 1 − XA ( τ)
1 + kτ
XA ( τ) =
(3.75)
Сопоставление формул (3.72)-(3.75) показывает, что при
равных значениях ХА величина Se в аппарате идеального
вытеснения окажется заметно выше, чем в аппарате идеального
смешения (рис. 3.24).
107
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Рис. 3.24. Зависимость селективности от степени
превращения реагента: 1 – идеальное вытеснение; 2 – идеальное
смешение
Если бы целевым компонентом было вещество С, то и в
этом случае селективность была бы выше для потока идеального
вытеснения.
Для параллельной реакции (В – целевой продукт)
k1
A 
→B
k2
A 
→C
влияние структуры потока на селективность зависит от
соотношения порядков основной и побочной стадии.
Если обе стадии одинакового порядка, то соотношение
скоростей образования обоих продуктов зависит только от
соотношения констант скорости:
rB k1C nA1 k1
=
=
.
rC k 2 C nA2 k 2
В этом случае селективность не меняется по ходу
реакции (если не меняется температура) и не зависит от типа
потока.
Теперь пусть первая стадия – первого порядка, а вторая –
второго. Рассчитаем соотношение rB к rC :
rB k1C A k1 1
=
=
.
rC k 2 C2A k 2 C A
Это отношение, а стало быть, и селективность тем
больше, чем меньше СА. Вспоминая рис. 3.23, видим, что в
аппарате идеального смешения концентрация компонента А
108
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
будет ниже, чем в аппарате идеального вытеснения. Значит,
процесс в аппарате идеального смешения проходит, хотя и
медленно, но зато с существенно большей селективностью.
Нетрудно понять, что если целевая реакция − второго порядка, а
побочная первого, то все преимущества будут на стороне
аппарата идеального вытеснения.
3.2.6. Влияние переноса массы на протекание
химической реакции и селективность
Важную роль в макрокинетике химических реакций
играют массообменные, или диффузионные факторы. Влияние
массообмена невелико, если реакция гомогенная: в этом случае
процесс идет в основной массе потока, где конвекция быстро
доставляет необходимые для реакции вещества.
По-другому обстоит дело в случае гетерогенных
процессов. Гетерогенные процессы проходят либо на
поверхности (например, гетерогенный катализ, кристаллизация,
электролиз, обжиг, горение твердого топлива), либо сопряжены
с переходом вещества из одной текучей фазы в другую
(абсорбция, экстракция, процессы в газожидкостных реакторах
и т.д.). Здесь реакция не может пройти без того, чтобы вещество
было доставлено к поверхности раздела фаз или через эту
поверхность из одной фазы в другую, поэтому непременная
стадия таких процессов — диффузия у границы раздела фаз.
Диффузия происходит гораздо медленнее, чем конвекция,
поэтому наличие диффузионной стадии может тормозить
химическую реакцию.
Закономерности
диффузионного
торможения
химических реакций сложны и многообразны. Многообразие
связано с особенностями подвода вещества к зоне реакции. Если
реакция протекает на наружной поверхности, омываемой
потоком, говорят о внешнедиффузионном торможении. Влияние
диффузии
при
реакции
в
порах
обусловливает
внутридиффузионное торможение. Здесь мы рассмотрим
закономерности внешнедиффузионного торможения на примере
каталитических реакций. Закономерности некаталитических
гетерогенных реакций в основном не отличаются от
рассматриваемых.
109
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Используя подход Д.Л. Франк-Каменецкого, будем
рассматривать простой случай: процесс стационарен, проходит
изотермически. Не углубляясь в механизм каталитической
кинетики, будем считать, что реакция необратима, имеет первый
порядок, а ее скорость равна
r = kC ПA ,
(3.76)
где
−
концентрация
реагирующего
вещества
CПA
непосредственно у поверхности.
Для реакций не первого порядка закономерности близки
к тем, которые получаются в анализируемом случае. Пусть
реакция протекает на непористом катализаторе. Поверхность
катализатора
омывается
потоком
газа;
концентрация
реагирующего вещества в ядре потока составляет CАЯ . В
результате протекания реакции концентрация у поверхности
уменьшается и составляет CПА . Разность (CАЯ − CАП ) является
движущей силой массоотдачи, в результате чего возникает
поток массы через диффузионный пограничный слой.
Рассмотрим материальный баланс по реагирующему
веществу. Скорость расходования реагента за счет реакции
составляет
 dM A 
П
rA = 
 = − kC А ,
 Fdt 
(3.77)
скорость подвода вещества диффузией равна
 dM A
jА = 
 Fdt

Я
П
 = β(C А − C А ) ,

(3.78)
где β — коэффициент массоотдачи; F – поверхность
массообмена.
Так как процесс стационарный, сумма прихода и расхода
равна нулю:
β(CАЯ − CАП ) = kC АП .
(3.79)
Откуда
110
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
CПА =
β C АЯ
.
β+k
(3.80)
Подставляя выражение (3.79) в формулу (3.77), получим
r = −rA =
kβ Я
CА = k 'C АЯ .
k +β
(3.81)
Уравнение (3.81) аналогично уравнению (3.78), но
вместо величины CП , которую мы не можем непосредственно
определить, в него входит величина C Я − концентрация в ядре
потока, а вместо истинной константы скорости реакции k −
эффективная константа k ' , определяемая по формуле (3.81).
Таким образом, эффективная константа скорости зависит
и от скорости реакции, и от скорости диффузии. Обратная k '
величина, называется общим сопротивлением процессу:
1 1 1
= +
k′ k β
(3.82)
Общее сопротивление равно сумме диффузионного и
кинетического сопротивлений. Из формулы (3.82) ясно, что
всегда k ' < k , то есть наличие диффузионного сопротивления
тормозит процесс.
Выделим некоторые характерные случаи.
1. Величина β много больше k . Это означает, что либо
очень велика скорость массоотдачи, либо очень мала скорость
реакции. В этом случае k ' ≈ k . Скорость процесса лимитирует
кинетика. Процесс происходит в кинетической области. При
этом уравнение (3.80) получает вид
CАЯ ≈ C АП ,
(3.83)
то есть диффузия успевает подводить к поверхности столько
вещества, что его концентрация у поверхности практически не
отличается от концентрации в ядре потока. Скорость реакции
равна
111
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
r ≈ kC АЯ .
(3.84)
2. Наоборот,
k
β . Скорость диффузии намного
меньше скорости реакции k ' β . Лимитирует массоотдача,
процесс происходит в диффузионной области. Тогда из
уравнения (3.80) следует
CПА ≈ 0 .
Все вещество, подводимое диффузией, успевает
полностью прореагировать в быстрой реакции. Для скорости
реакции получается
r ≈ β CАЯ .
(3.85)
3. Коэффициенты k и β сравнимы по величине. Ни одна
из стадий не лимитирует процесса - он происходит в
промежуточной (переходной) области. Скорость реакции
необходимо рассчитывать по уравнению (3.81).
Увеличению скорости массоотдачи в первую очередь
способствует рост скорости потока. В то же время скорость
потока непосредственно не влияет на химическую кинетику,
поэтому, если процесс происходит в кинетической области,
гидродинамика не влияет на его скорость. Если в
диффузионной, то увеличение скорости потока ускоряет
процесс. Если в промежуточной, то при увеличении скорости
потока процесс может перейти в кинетическую область.
На скорость химической реакции сильнее всего
оказывает влияние температура. В то же время скорость
массоотдачи слабо зависит от температуры, в особенности если
движущая фаза − газ. В связи с этим кривая зависимости
скорости реакции от температуры имеет характерный вид (рис.
3.25).
В области низких температур общая скорость процесса
сравнительно мала, но при этом с повышением температуры она
резко растет по экспоненте в соответствии с уравнением
Аррениуса. Процесс идет в кинетической области. При
повышении температуры рост скорости процесса замедляется и,
начиная с некоторой температуры, почти прекращается.
Процесс вступает в диффузионную область.
112
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Рис. 3.25. Зависимость скорости реакции от температуры
Если увеличить скорость потока, то есть снизить
диффузионное сопротивление, то в кинетической области
кривая пойдет точно так же, но промежуточная и диффузионная
области начнутся при несколько более высоких температурах, а
общая скорость процесса в этих областях окажется выше.
Как уже говорилось, описанные закономерности верны
не только для катализа на непористой поверхности, но и для
многих других гетерогенных процессов. Кинетика ряда
процессов
горения
и
обжига
определяется
внешнедиффузионным торможением. С ним связаны явления
концентрационной поляризации в электролизе.
Влияние диффузионной кинетики на селективность
иногда очень существенно и может определять успех или
неуспех того или иного варианта технологии.
Селективность параллельной реакции
1 по А; k1
A 
→В
1 по А; k 2
A 
→С
зависит лишь от отношения констант скоростей. Диффузионное
торможение не влияет на селективность.
Иное дело – последовательная реакция (целевой продукт
– В)
1 по А; k1
A 
→B
Аналогично
зависимости
1 по В; k 2
B 
→C
уравнению (3.79),
113
легко
получить
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
βA (CAЯ − CAП ) = k1CAП

βB (CBП − CBЯ ) = k1C AП − k 2 C BП
β (C − C ) = k C
CЯ
2 BП
 C CП
(3.86)
В уравнениях (3.86) отражено направление переноса
вещества: к поверхности – для А, в ядро потока – для В и С.
Отсюда легко получить
β C
CAП = A AЯ ;
β A + k1
CBП =
βA C AЯ k1 + ( k1βB + βAβB ) C BЯ
.
( k1 + βA )( k 2 + βB )
В самом начале процесса CВЯ =0, тогда выражение для
CВП примет вид:
CBП =
Тогда получим
βA CAЯ k1
.
k
+
( 1 βA )( k 2 + βB )
rB β B
.
=
rC k 2
(3.87)
Из формулы (3.87) следует, что переход второй стадии
реакции из кинетической области в диффузионную коренным
образом изменит ход реакции. В кинетической области
(βB / k 2 ) → ∞ , и в начале процесса Se=1. Образующийся
продукт В столь быстро отводится от поверхности, что вторая
стадия вообще вначале не идет, начинаясь лишь тогда, когда В
накопится в ядре потока. В диффузионной же области
(βB / k 2 ) ≈ 0 и S=0, причем селективность нулевая на всем
протяжении процесса: целевой продукт не успевает
диффундировать в поток и практически целиком вступает в
побочную реакцию. В этом случае совершенно необходимо
интенсифицировать
массообмен,
снять
диффузионное
торможение и работать в кинетической области.
114
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
IV. Методы и модели определения физикохимических свойств газовых и жидких смесей
4.1. Уравнения состояния
Для
моделирования
промышленных
процессов
химической технологии необходима надежная информация по
различным физико-химическим свойствам индивидуальных
веществ и их смесей. Эта информация составляет обычно более
половины требуемых исходных данных, и от ее надежности во
многом зависит точность и адекватность получаемых
результатов. В виду того что экспериментальная база данных, по
понятным причинам ограничена, в инженерных расчетах для
определения физико-химических свойств рабочих агентов
используют модели. Эти модели, как правило, в своей основе
имеют либо корреляционные зависимости, построенные
посредсвом обобщения большого числа экспериментальных
данных, либо некоторую полуэмпирическую теорию,
базирующуюся на упрощенном представлении физической
картины рассматриваемого явления, поэтому такие модели
имеют ограничения по области применения и нуждаются в
экспериментальном определении соответствующих параметров.
Кроме
того,
основным
недостатком
известных
полуэмпирических моделей является то, что каждая из них
способна
описывать
только
одну
конкретную
термодинамическую характеристику, для определенного класса
веществ или смесей, поэтому задача построения надежного,
универсального метода расчета физико-химических свойств
газов и жидкостей является и в настоящее время актуальной.
Идеальный метод расчета физико-химических свойств
газов
и
жидкостей
должен
обладать
следующими
характеристиками:
1. Выдавать надежные данные для чистых веществ и их
смесей во всей области термодинамических состояний, включая
фазовые границы.
2. Указывать агрегатные состояния вещества.
3.
Требовать
минимальное
количество
экспериментальных данных.
4. Обладать невысокой вычислительной трудоемкостью.
115
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Все физико-химические свойства можно разделить на
равновесные (давление P, внутренняя энергия U, удельная
теплоемкость Cp, r – теплота испарения, плотность ρ) и
неравновесные (коэффициент молекулярной теплопроводности
λ, коэффициент кинематической молекулярной вязкости ν,
коэффициент молекулярной диффузии D). В разделе 3.2 было
показано, что расчет равновесных свойств веществ не
представлял бы никаких проблем, если было бы известно
уравнение состояния
P = f (T, v, x1 , x 2 ,K x n −1 ) . Строгую
математическую формулировку уравнения состояния удалось
получить только в рамках молекулярно-статистической теории в
середине прошлого столетия. Однако в настоящее время ее
использование во многих случаях затруднительно ввиду
сложностей математического характера.
Исторически уравнение состояния пытались получить,
обобщая экспериментальные данные по поведению газов и
жидкостей, а в дальнейшем привлекая теорию молекулярного
строения веществ.
Первым уравнением состояния является уравнение
идеального газа Менделеева−Клапейрона, построенное на
основе обобщения экспериментальных результатов Бойля,
Клайперона, Шарля, и др.:
PV = NRT .
Это уравнение оказалось справедливым только для газов
умеренной плотности, то есть такого состояния, в котором
молекулы вещества не «чувствуют» друг друга, поэтому оно с
большой ошибкой описывало поведение плотных газов и не
предсказывало
существование
фазовых
переходов.
Экспериментальные данные показывали, что диаграмма
состояния веществ имеет характерный вид, представленный на
рис.4.1.
На ней присутствуют линии конденсации (ОК), кипения
(КВ) и кристаллизации (ВМ), а также характерные точки:
критическая К и тройная TP. Существенным шагом вперед стало
уравнение, предложенное Ван-дер-Ваальсом (ВдВ) в 1908 г. и
носящее его имя:
116
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
P=
RT
a
− 2,
VM − b VM
(4.1)
где параметр b учитывает собственный размер молекул,
параметр a – притяжение молекул друг к другу, в результате
которого давление в системе уменьшается; VM − молярный
объем.
Рис. 4.1. Диаграмма состояния
117
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Рис. 4.2. Изотермы Ван-дер-Ваальса
На ней присутствуют линии конденсации (ОК), кипения (КВ) и
кристаллизации (ВМ), а также характерные точки: критическая
К и тройная TP. Существенным шагом вперед стало уравнение,
предложенное Ван-дер-Ваальсом (ВдВ) в 1908 г. и носящее его
имя:
P=
RT
a
− 2,
VM − b VM
(4.2)
где параметр b учитывает собственный размер молекул,
параметр a – притяжение молекул друг к другу, в результате
которого давление в системе уменьшается; VM − молярный
объем. Уравнение ВдВ − это кубическое уравнение по объему,
которое при температурах, ниже критической, имеет 3 корня,
два из которых (наименьший и наибольший) соответствуют
объемам равновесных паровой и жидкой фаз (рис. 4.2). При
температуре, выше критической, уравнение ВдВ имеет один
корень.
Таким
образом,
уравнение
ВдВ
позволяет
118
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
количественно определять существование жидкой и паровой
фаз, а также критическое состояние вещества. Значения
констант в уравнении (4.2) можно связать с критическими
параметрами
вещества
и
Tc ,
Pc
vc , используя
термодинамические условия поведения давления в критической
точке:
∂ 2 P(T, VM )
∂P(T, VM )
=0 и
=0,
∂VM
∂VM 2
откуда
vc = 3b Pc =
a
8a
, Tc =
.
2
27b
27Rb
Как оказалось, точности уравнения ВдВ не достаточно
для описания плотных газов и жидкостей. Тем не менее, удачная
структура уравнения ВдВ была положена в основу большого
количества уравнений состояния, широко используемых в наши
дни, например уравнения Редлиха−Квонга:
P=
RT
a
+ 0.5
.
VM − b T VM (VM + b)
Важным следствием уравнения ВдВ стала теория
соответственных состояний, согласно которой в приведенных
переменных все вещества подчиняются одному и тому же
уравнению состояния. За приведенные переменные обычно
принимают
P* =
T
P
V
, T* =
, v* = M .
Tc
Pc
vc
Уравнение ВдВ, выраженное в этих переменных,
принимает вид
P* =
8T *
3
− 2.
3v * −1 v *
Теория соответственных состояний подтверждается для
большого класса веществ, даже если уравнение для них не
совпадает по форме с (4.2), однако она не оказалась общей, и
попытки найти универсальное уравнение состояния в
приведенных переменных не увенчались успехом.
119
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Важное место среди моделей уравнений состояния
занимает вириальное уравнение
P=
RT B2 (T) B3 (T)
+
+
+K
VM
VM 2
VM 3
(4.3)
Первоначально это уравнение было предложено Тизеном
(1885 г.) как чисто эмпирическое и долгое время было
предметом детального исследования многих ученых. Позднее
Урсел (1927 г.) показал, что у этого уравнения существуют
строгие молекулярно-статистические основы. Смысл уравнения
(4.3) заключается в том, что давление может быть представлено
в виде ряда по плотности с коэффициентами, которые зависят
только от температуры. Впоследствии было доказано, что для
жидкофазного состоянии, где плотность высока, этот ряд
расходится. Тем не менее, существует большое количество
модификаций уравнения (4.3), которые используются в
практических расчетах, среди которых наиболее известное
уравнение Бенедикта−Вебба−Рубина:
C0 2
)ρ −
T2
.
cρ3
3
6
2
− γρ2
−(bRT − a)ρ + aαρ + 2 (1 + γρ )e
T
P = RTρ + (B0 RT − A 0 −
(4.4)
Уравнение (4.4) может быть расширено в зависимости от
желаемой расчетной области. Например, для индивидуальных
предельных углеводородов уравнение типа (4.4), способное
описывать области пара, жидкости и небольшую область
сверхкритических состояний, имеет 32 настраиваемых
параметра [13].
Таким образом, можно констатировать, что в настоящее
время не существует единого уравнения состояния, поэтому при
выборе того или иного варианта нужно обратить особое
внимание на область применимости данного уравнения. Даже
уравнения, содержащие 30−40 констант, имеют значительные
ограничения, а попытки обобщить уравнения с тем, чтобы они
охватывали более широкий круг веществ, обычно приводят к
снижению точности. Сравнение большого числа уравнений
120
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
состояния и предварительные рекомендации по их выбору
представлены в [14].
В заключение этого раздела обратим внимание на
очевидный факт: в химической технологии рабочими агентами
являются многокомпонентные смеси. В этом случае в уравнении
состояния возникают дополнительные переменные, связанные с
составом: p = f (T, v, x1 , x 2 ,K x n −1 ) . Очевидно, что получить
обобщение представленных выше уравнений состояния на
смеси можно только в том случае, если компоненты смеси
описываются
однотипными
уравнениями
состояния.
Практически такое обобщение осуществляется посредством
выражения параметров уравнения состояния в виде функции от
состава, например для уравнения ВдВ необходимо записать
и
a m = f1 (x1 ,K , x n −1 )
b m = f 2 (x1 ,K , x n −1 ) .
(4.5)
Соотношения типа (4.5) называются правилами
смешения, и от их вида очень сильно зависит поведение
макросвойств смеси в зависимости от состава. Теоретически
правила смешения получить невозможно, поэтому все они
являются эмпирическими и сформулированы в результате
многочисленных проб и сравнений с экспериментом. Наиболее
часто используются следующие соотношения:
n
n
n
a m = ∑ x i a ik , a m = ∑∑ x j x i a ijl ,
i =1
i =1 j=1
где аi и аij – унарные и бинарные параметры. Иногда используют
выражения и с тройными параметрами. Для определения
бинарных параметров обычно применяют следующее правило:
a ij = ηij
ai + a j
2
− если параметр имеет физический смысл
геометрического размера;
a ij = ξij a i a j − для параметров, отражающих энергетическую
составляющую.
4.2. Расчет термодинамических свойств на основе
избыточных функций
121
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Этот метод используется в том случае, когда расчет на
основе уравнения состояния невозможен или не дает требуемой
точности результатов.
Функцией
смешения
Am
называют
изменение
термодинамической функции А при образовании раствора из
чистых компонентов:
n
A m (P, T, x1 , x 2 ,...x n −1 ) = A(P,T, x1 , x 2 ,..., x n −1 ) − ∑ x i Ai0 (T, P).
i =1
В частности, для энергии Гиббса:
x i µi − ∑ x i µi0
Gm
∑
=
= ∑ x i ln(x iγ i ) .
NRT
RT
Избыточной термодинамической функцией AE называют
разность между функций смешения реального и идеального
растворов:
A E (P, T, x1 , x 2 ,...x n −1 ) = A m (P, T, x1 , x 2 ,...x n −1 ) −
−Aidm (P, T, x1 , x 2 ,...x n −1 )
.
Для энергии Гиббса
GE
= ∑ x i ln(x iγ i ) − ∑ x i ln(x i ) =∑ x i ln(γ i ) .
NRT
Таким образом, видно, что избыточную энергию Гиббса
составляют важные для расчета условий равновесия
коэффициенты активности, поэтому вместо определения их
через уравнение состояния можно попытаться построить модель
непосредственно для избыточной энергии. Коэффициенты
активности
в
этом
случае
определятся
согласно
термодинамическому соотношению
1 ∂G E
ln( γ i ) =
RT ∂N i
.
p,T,Ni
Такой подход может применяться для определения не
только коэффициентов активности, но и любой другой
термодинамической функции. Недостатком в этом случае будет
122
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
отсутствие термодинамической связи между моделями для
различных величин.
Для построения модели избыточной энергии Гиббса
нужно использовать два строгих соотношения:
- асимптотического поведения G E = 0 при x1 → 0 и
x2 → 0 ;
- Гиббса−Дюгема (3.5).
Рис. 4.3. Различные зависимости избыточной функции от
концентрации бинарной системы
Таким образом, вид кривой избыточной функции может
иметь самый различный характер, но в крайних точках (для
чистых веществ) она должна быть равна нулю, как это следует
из самого понятия избыточной функции (рис. 4.3), а получаемые
из избыточной энергии Гиббса коэффициенты активности
должны удовлетворять соотношению Гиббса−Дюгема.
Одна из самых ранних и простых моделей для
избыточной энергии Гиббса – модель Маргулиса:
GE
= A12 x1x 2
NRT
и
123
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ln( γ1 ) = A12 x 22
ln( γ 2 ) = A12 x12
,
здесь параметр бинарного взаимодействия A12 определяется из
эксперимента.
Эта модель достаточно точно описывает поведение
систем с небольшим отклонением от идеальности. Также
простые аппроксимационные формулы для коэффициентов
активности могут иметь вид степенных рядов:
n
ak k
x2
k =2 k
.
m
bk k
ln( γ 2 ) = ∑ x1
k =2 k
ln( γ1 ) = ∑
На практике часто используются модели, полученные на
основе упрощенного рассмотрения физической картины
явлений, происходящих в смеси. Наиболее распространены
модели:
• решеточные и ячеечные;
• регулярных растворов;
• локального состава;
• групповых составляющих
и их различные комбинации.
Широкое распространение получили уравнения теории
локального состава, которые предназначены прежде всего для
расчета коэффициентов активности компонентов жидких
растворов. Концепция локального состава впервые была
сформулирована в работе Вильсона в виде соотношений,
приближенно описывающих локальную упорядоченность
молекул в растворе. В результате была предложено следующее
уравнение Вильсона:
k
 k

GE
= −∑ x i ln  ∑ x i λ ij 
NRT
i =1
 j=1

λ ij =
Vj
Vi
e
−
Cij
RT
(4.6)
124
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»






x λ
ln( γ i ) = 1 − ln  ∑ x jλ ij  − ∑  k m mi 
 j=1
 m =1  ∑ x λ 
j mj 

 j=1

k
k
где Cij − параметр, характеризующий энергию взаимодействия
молекул сортов i и j , Vi ,VJ − мольные объемы компонентов.
Параметры
этого
уравнения
находятся
по
экспериментальным данным о равновесии пар−жидкость в
бинарной смеси. Ясно, что при этом параметры оказываются
привязанными к конкретным веществам. Уравнение Вильсона
не пригодно для описания равновесия жидкость−жидкость.
Развитием данной теории стали модель NRTL, UNIQUAC,
возможности которых были расширены на описание равновесий
жидкость-жидкость. Хороший обзор наиболее известных
моделей для расчета коэффициентов активности, а также
рекомендации по их выбору представлен в [4,14]. При
существующем большом многообразии моделей все они имеют
ограниченную область применения, как по физическим
системам, так и по области состояний, поэтому необходимо
проявлять некоторую осторожность при выборе той или иной
модели. Кроме того, основной задачей большинства моделей
является не предсказание термодинамических свойств, а лишь
их аппроксимация. Единственное, что часто удается, это
описать свойства многокомпонентных систем по известным
свойствам бинарных, см., например, (4.6), поэтому для
использования этих моделей обязательно нужны надежные
экспериментальные данные, по крайней мере, по бинарным
системам. Учитывая ограниченность экспериментальной базы
данных, что естественно в условиях, когда в год синтезируются
тысячи новых соединений, важно иметь модели способные
предсказывать поведение термодинамических свойств чистых
веществ и смесей по их молекулярному строению.
Подобная идея реализована в моделях групповых
составляющих. Здесь рассчитываемое свойство представляется
как сумма вкладов атомов и атомных группировок, из которых
состоит молекула. Если принять, что вклад одних и тех же
125
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
атомов в разных соединениях будет одинаков, то, учитывая, что
число атомов и атомных группировок ограничено, можно
получить некоторый небольшой унифицированный набор
параметров, с помощью которого можно будет рассчитывать
свойства различных веществ. К сожалению, гипотеза об
одинаковом вкладе атомов в различных молекулах оказалась
несостоятельной, и для большей точности расчета нужно
учитывать влияние соседних атомов и группировок. Однако
даже в этом случае число параметров, необходимых для расчета
свойств оказывается несоизмеримо меньше того количества
веществ и их смесей, которые можно на их основе рассмотреть.
Метод групповых составляющих был реализован в модели
UNIFAC и различных ее модификациях. Эта модель является
основной в инженерных методах расчета и в системах
автоматизированного
проектирования,
когда
экспериментальные данные по одному или нескольким
компонентам рассматриваемой системы отсутствуют.
4.3. Молекулярно-статистические методы описания
макросвойств газов и жидкостей
Известно, что наблюдаемые в физическом эксперименте
макросвойства веществ связаны с их молекулярным строением.
Молекулярное строение определяет характер и величину силы
взаимодействия между двумя отдельными молекул, что в свою
очередь определяет коллективное поведение большого числа
подобных
молекул
и
соответственно
наблюдаемые
макросвойства вещества. Основной задачей молекулярностатистической теории является установление связи между
взаимодействиями отдельных частиц, рассматриваемой системы
и макросвойствами, характеризующими систему как целое.
Существование такой связи в виде математических
соотношений, решение которых на современных ЭВМ занимало
бы разумное время, являлось бы наиболее эффективным и
перспективным методом описания свойств газов и жидкостей.
Исходными данными в этом случае являются потенциалы
межмолекулярного
взаимодействия
компонентов
рассматриваемой системы. Потенциал межмолекулярного
взаимодействия зависит от атомного строения молекулы
126
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
вещества и незначительно − от параметров состояния, что
позволяет говорить о потенциалах межмолекулярного
взаимодействия как об универсальной базе данных. Причем эти
потенциалы, как и геометрия молекулы, могут быть рассчитаны
методами квантовой химии.
Размеры молекул реальных веществ оказываются
достаточно большими для того, чтобы в хорошем приближении
их динамическое поведение можно было описывать на основе
законов классической механики. Тогда состояние системы из N
одинаковых молекул можно определить заданием координат
r
r
{r} и импульсов {p} всех молекул, которые в данном случае
являются каноническими переменными. Мгновенные значения
канонических переменных определяют микросостояние системы
и ее микросвойства. Изменение координат и импульсов
отдельных молекул системы подчиняется уравнениям
движения:
r ∂H
ri = r
∂ pi
r
∂H
pi = −
,
∂ri
,
i=1, 2,K , N
(4.7)
где H − гамильтониан систем, которая определяется как сумма
кинетической и потенциальной энергии молекул:
N r2
rr
r
p
H(r, p) = ∑ i +U(r) ,
i =1 2m
(4.8)
где m − масса частицы. Потенциальную энергию можно
представить в виде ряда
U N (r1 , r2 ,K , rN ) = ϕi (ri ) + ϕij (ri , rj ) +
ϕijk (ri , rj , rk ) + K ,
∑
(4.9)
где ϕi − потенциал
∑
i< j
внешнего
∑
i < j< k
поля;
ϕij
−
потенциал
межмолекулярного взаимодействия. Для большинства реальных
систем достаточно учесть только бинарные потенциалы
межмолекулярного взаимодействия.
127
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Все измеряемые прибором макросвойства веществ
являются средними по времени измерения, величинами от
динамических функций системы. Под динамическими
функциями понимают функции, определяющие мгновенные
значения микросвойств системы в зависимости от канонических
переменных. Время изменения динамических переменных
намного меньше времени измерения прибором, поэтому
измеряемые макропараметры можно представить как
D
τ
τ
r
r
1
= ∫ D(r(τ), p(τ))dτ ,
τ0
(4.10)
r
r
где D(r(τ), p(τ)) − динамическая функция системы,
соответствующая измеряемому макросвойству; τ − время
наблюдения.
Таким образом, решая уравнения (4.7) можно рассчитать
траектории движения молекул, составляющих рассматриваемую
систему, и по известной динамической функции определить
макросвойство системы по (4.10). Такой способ расчета
макросвойств веществ называется молекулярной динамикой
(МД). Очевидно, что подобная задача может быть решена лишь
с помощью ЭВМ, причем, для того чтобы макросвойства
рассматриваемой в численном эксперименте системы были
близки к реальной, число молекул N должно быть достаточно
велико. В идеале значение N должно иметь порядок 1023,
соответствующий числу молекул в 1 моле вещества.
Современные ЭВМ позволяют рассмотреть систему, состоящую
из N 103− 4 . В принципе этого оказывается достаточно, так как
значения макросвойств, близкие к реальным, начинают
проявляться у систем состоящих из нескольких сотен молекул.
Каждому
макроствойству
можно
поставить
в
соответствие свою динамическую функцию. Например, для
определения макроскопической плотности в точке x выделим
некоторую малую область пространства G в ее окрестности.
Введем функцию f G (r) таким образом, что:
f G (ri ) = 1 , ri ∈ G ;
f G (ri ) = 0 , ri ∉ G .
128
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
В пределе точечных частиц будем иметь δ функцию
Дирака f G (ri ) = δ(x − ri ) . Тогда микроплотность частиц в точке
x будет равна
n(x, r) =
1 N
∑ fG (ri ) ,
G i =1
(4.11)
а макроплотность получится усреднением соответствующей
динамической функции (4.11):
ρ(x) = n(x, r(τ)) τ .
Динамической функцией для выражения внутренней
энергии системы является сам гамильтониан (4.8), тогда
r
r
U = H(p(τ), r(τ)) .
τ
(4.12)
Для определения давления удобно воспользоваться
теоремой о вириале сил. Пусть рассматриваемая система,
состоящая из N молекул, заключена в объеме прямоугольной
формы при температуре T (рис. 4.4) . Согласно теореме о
вириале сил
N
∑rf
i =1
i i
=2 K
(4.13)
129
, K =
3
Nk B T .
2
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Рис. 4.4. К теореме о вириале сил
Разделим силы, действующие на молекулы системы, на
внутренние fv, создаваемые другими молекулами, и внешние fs,
которые являются следствием действия на систему
ограничивающих ее стенок. Внешние силы, очевидно, можно
представить через величину внешнего давления:
N
∑ r fs
i =1
i
= −(L z PSxy + L x PSyz + L y PSxz ) = −3PV .
i
(4.14)
Внутренние
силы
запишем
межмолекулярного взаимодействия:
N
N
N
N
N
∑ r fv = ∑ r ∑ fv = ∑∑ r fv
i =1
i
i
i
i
ij
j
i
i< j
ij
ij
через
потенциал
N
N
∂ϕ (rij )
i
i< j
∂rij
= − ∑∑ rij
.
(4.15)
Здесь было учтено, что согласно закону Ньютона
fv ij = −fv ji ,
а также rij = rj − ri − межмолекулярное расстояние. Подставляя
(4.14), (4.15) в (4.13), получим
pV = Nk B T −
1
3
N
N
∑∑ r
i
i< j
ij
∂ϕ (rij )
∂rij
.
(4.16)
Общая
расчетная
схема
термодинамических
характеристик методом МД выглядит следующим образом:
1. Задается состояние системы: число частиц,
температура, плотность или давление, а также вид потенциала
межмолекулярного взаимодействия ϕij .
2. Фиксируется начальное микросостояние, то есть
начальные значения координат и импульсов r0 и p 0 всех
частиц. Обычно за начальное выбирается упорядоченное
расположение молекул в узлах кристаллической решетки.
Начальные скорости можно задать случайным образом в
соответствии с условием, чтобы суммарный импульс системы
был равен нулю и значением температуры.
130
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
3. На определенном временном интервале решается
система 6N уравнений (4.7). Для этого она представляется в
конечных разностях.
4. По рассчитанным траекториям движении молекул
определяются различные макросвойства системы, например
внутренняя энергия по(4.12), давление по (4.16) и т.д.
Ошибку, вносимую в расчет конечностью числа частиц,
частично
компенсируют
периодическими
граничными
условиями. Если в процессе моделирования частица выходит
через какую-либо грань ячейки, то точно такая же частица
заходит через противоположную грань. Таким образом, удается
сохранить число частиц и энергию системы.
Метод МД имеет большие возможности, так как
позволяет описать изменение наблюдаемых величин, как в
пространстве, так и во времени, что дает возможность находить
не только равновесные, но и кинетические свойства системы.
Метод ансамблей Гиббса
Вместе с детерминированным подходом описания
движения молекул уравнениями Ньютона существует модель
молекулярного хаоса, сформулированная Больцманом для
описания теплового движения молекул. Согласно этой
стохастической модели движение молекул из-за свой сложности
носит случайный (вероятностный) характер. Именно эта модель,
как правило, ложится в основу качественного анализа явлений,
протекающих за счет теплового движения молекул, например
диффузии, теплопроводности и т.д. Каким же образом
согласуются эти два подхода: возможность точного расчета
положения и импульса молекул в любой момент времени и
утверждение о молекулярном хаотическом движении? Если
проанализировать сущность метода МД, то станет ясно, что
единственным местом в этой теории, где может возникнуть
случайная компонента, являются начальные условия. Системы,
состоящие из большого числа молекул, являются примером
неустойчивых систем, для которых небольшое отклонение в
начальных условиях может приводить к значительному
изменению рассчитываемых траекторий. Попробуем оценить
требуемую точность задания начальных условий, если мы хотим
131
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
проследить за траекторией одной частицы в течение одной
секунды. Будем считать эту попытку успешной, если ошибка в
определении направления движения не будет превышать угол π,
то есть мы сможем всего лишь определить, летит частица
вперед или назад. Предположим, что две сферические молекулы
газа сталкиваются под углом α, причем этот угол определен с
точность δ 0 . Будем считать столкновения молекул абсолютно
упругими. Ввиду того что поверхность сталкивающихся частиц
является выпуклой, ошибка в определении угла после каждого
столкновения будет увеличиваться, и после k столкновений ее
можно оценить как
δ 
δk = δ0  1 
 δ0 
k
или, определяя k = t / τ , где t – время; τ − среднее время между
двумя
 t
δ(t) = δ0 exp  α  ,
 τ
столкновениями,
где
δ 
α = ln  1  ≈ 1 . Поскольку для газов при нормальных условиях
 δ0 
τ 10 −12 , для выполнения поставленных условий начальный
угол должен быть определен с точностью δ0 = exp ( −1012 ) , что
далеко выходит за все разумные пределы. Таким образом,
несмотря на то, что в молекулярных системах частицы движутся
по законам классической механики, траектории их движения
можно считать непредсказуемыми.
Основываясь на вероятностном представлении теплового
движения молекул, Гиббсом был предложен метод построения
описания макросвойств веществ. Здесь наблюдение за
поведением одной системы во времени заменяется наблюдением
за состояниями большого количества идентичных систем
(ансамбля) в некоторый момент времени. Пусть имеется
некоторое число одинаковых систем (число молекул, величина и
форма объема, температура полностью совпадают). Такая
совокупность систем называется ансамблем Гиббса. Задавшись
в каждой системе начальными условиями, зафиксируем
микросостояние каждой системы через определенный
132
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
промежуток времени. Эти состояния будут отличаться друг от
друга, так как задать идентичные начальные условия вследствие
неизбежной ошибки невозможно. Если ввести понятие 6Nмерного
фазового
r
r
r
r
r
r
объема: dГ = dr1 , dr 2 ,K , dr N , dp1 , dp 2 ,K , dp N ,
то
системе будет соответствовать своя точка в
пространстве (рис. 4.5).
каждой
фазовом
Рис. 4.5. Ансамбль Гиббса
При N → ∞ фазовый объем полностью заполнится
точками состояний. Тогда вероятность нахождения системы в
некоторой элементарной части фазового объема будет
пропорциональна:
W(r, p) = F (r, p)drdp ,
где F (r, p) − функция распределения или плотность
вероятностей.
Если известен вид функции распределения, то средние
по ансамблю значения динамической функции системы
определяются следующим образом:
D
(4.17)
q
= ∫ D(r, p, τ)F (r, p, τ)drdp .
Γ
Важнейшим постулатом статистической физики является
утверждение, что средние по времени значения (4.10)
динамических функций совпадают со средними по ансамблю
значениями (4.17) этих же динамических функций. Обоснование
133
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
этой гипотезы является предметом эргодической теории, и в
настоящее время не существует строгого ее доказательства. Из
основного постулата следует, что для нахождения наблюдаемых
на практике величин, а также их производных по времени и
пространственным координатам, необходимо найти явный вид
функции распределения.
Вид функции распределения зависит от типа системы
(закрытая, изолированная и т.д.), характера переменных,
определяющих состояние системы. Подробно ознакомиться с
различными ансамблями и соответствующим им функциями
распределения можно в монографии [15]. Например, для так
называемого канонического ансамбля, когда задаются N, V, T ,
F (r, P) = exp(− H / kT + F / kT)
и
W(r, P) = exp(− H / kT + F / kT)drdp ,
где А − некоторая константа, определяемая из условия
нормировки:
∫ exp(−H / kT + F / kT)drdP = 1 ;
(4.18)
F = − kT ln(Z N ) ,
(4.19)
где
Z N = ∫ exp(− H / kT)drdP
−
статистическая
сумма
(интеграл).
Нетрудно показать, что F является свободной энергией
Гельмгольца.
Продифференцируем (4.18) по T :
 H
∫  kT
2
+
∂  F 

 exp(− H / kT + F / kT)drdP = 0 ;
∂T  kT  
∂  F 
1
H exp(− H / kT + F / kT)drdP =

.
2 ∫
∂T  kT 
kT
Интеграл в левой части
термодинамическое выражение для
системы:
134
представляет собой
внутренней энергии
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
∂  F 
1
U=

.
2
∂T  kT 
kT
(4.20)
Это
известное
термодинамическое
соотношение
внутренней и свободной энергиями.
Таким образом, существует возможность определять
макросвойства системы не только посредством усреднения
соответствующих динамических функций, но и по вычисленной
статистической сумме
Z N . Действительно, вычислив
статистическую сумму, по (4.19) определяем свободную
энергию и далее, используя термодинамические соотношения
Максвелла, любые термодинамические свойства. Необходимо
заметить, что размерность задачи вычисления статистической
суммы не стала меньше, чем при усреднении динамических
функций по времени. Видно, что размерность интеграла
оказывается пропорциональной числу частиц в системе.
по импульсам выполняется
Интегрирование Z N
достаточно просто, учитывая вид гамильтониана (4.8) и
независимость друг от друга импульсов разных молекул,
получим
3
N
 mkT  2
ZN = 
 QN ,
 2πh 
(4.21)
∫
где Q N = exp(− U N (r1 , r2 ,K , rN ) / kT)dr − конфигурационный
интеграл.
Множитель
в
(4.21)
определяет
вклад
в
термодинамические свойства идеальногазовой составляющей.
Учитывая (4.20) и (3.9), определим внутреннюю энергию и
давление следующим образом:
U = kT 2
∂ (ZN )
∂ (Q N )
3
= NkT + kT
∂T V,N 2
∂T V,N
P = kT
∂ (Z N )
∂ (Q N )
= kT
∂V T,N
∂V T,N
135
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Основной
задачей
здесь
является
расчет
конфигурационного интеграла, для чего необходимо провести
интегрирование по координатам всех частиц системы. Учитывая
их количество в одном моле вещества, подобная задача
представляется невычислимой. Однако, как и в случае
молекулярно-динамического
моделирования,
остается
возможность рассмотреть систему, содержащую несколько
сотен или тысяч молекул. Конечно, даже в этом случае
использовать для интегрирования известные квадратурные
алгоритмы будет не эффективно. Наиболее приемлемым
методом решения таких задач является использование метода
Монте-Карло (МК). Основой идей этого метода является замена
вычислительных операций на упорядоченной расчетной сетке
вычислениями на случайным образом сгенерированных точках
расчетной области. При определенных условиях генерирования
этих случайных точек число вычислений для получения
конечного результата может оказаться намного меньше, чем в
квадратурных алгоритмах.
Для задач молекулярно-статистической теории метод МК
используется в варианте, предложенном Метрополисом [16].
Основой данного варианта является процедура, в которой
вместо случайного выбора конфигураций (микросостояния)
системы организуется цепь переходов системы по состояниям,
появление
которых
пропорционально
exp(−U N (r1 , r2 ,K , rN ) / kT ) . Алгоритм этого метода можно
представить следующим образом:
1. Задается состояние системы (число частиц,
температура, плотность или давление, а также вид потенциала
межмолекулярного взаимодействия ϕij ).
2. Фиксируется начальное микросостояние, т.е.
начальные значения координат r 0 всех частиц. Обычно за
начальное выбирается упорядоченное расположение молекул в
узлах кристаллической решетки. Определяется значение
потенциальной энергии системы в начальном состоянии. Если
учитываются только двухчастичные взаимодействия в
отсутствии внешних полей, то
136
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
U 0 (r10 , r20 ,K , rN0 ) =
1
ϕij (ri0 , rj0 ) .
∑
2 i< j
3. Генератор случайных чисел (ГСЧ) выдает три
случайных числа ξ1 , ξ 2 и ξ3 . По первому выбирается
молекула, которая будет двигаться, по второму и третьему –
направление ее движения (например, вдоль какой оси и в каком
направлении). Вычисляется новое положение выбранной
молекулы rξ11 = rξ01 + ∆r и соответствующее ему значение
энергии системы : U1 = U(r1 ) .
4. Проверяется условие приемлемости выполненного
перемещения. Если U1 < U 0 , то новое положение молекулы
принимается и алгоритм повторяется заново с п. 3. Если
U1 > U 0 , ГСЧ выдает следующее число ξ 4 , и проверяется новое
условие exp(−(U1 − U 0 )) < ξ 4 . При положительном результате
новое положение молекулы принимается, и алгоритм
повторяется с п. 3, при отрицательном − молекула возвращается
в исходное положение, алгоритм повторяется с п. 3.
Среднее значение динамической функции по такой
цепочке передвижений получают простым суммированием ее
значений в каждой конфигурации. Причем, если какой-либо
переход не был принят, значение динамической функции в
старом состоянии учитывается повторно:
U =
1 M
U(r i ) ,
∑
M i=1
где М − число сгенерированных конфигураций.
Для выполнения условия постоянства числа частиц, так
же как в методе молекулярной динамики, используют
периодические граничные условия.
Метод МК оказывается удобнее по сравнению с методом
молекулярной динамики для чисто изотермических процессов,
так как температура в этом методе легко фиксируется. Кроме
того, он также легко применяется как для различных типов
ансамблей (NVT), или (NPT) ансамблей с переменным числом
частиц, так и для систем, в потенциальной энергии которых
необходимо учесть не только парные межмолекулярные
137
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
взаимодействия. К недостаткам этого метода можно отнести
невозможность определения неравновесных свойств.
Методы МД и МК в силу их фундаментальной строгости,
а также надежности, получаемых на их основе результатов,
называют методами численного эксперимента. Очевидно, что
возможности этих методов напрямую зависят от современного
развития вычислительных средств, а также от наличия
реалистичных потенциалов межмолекулярного взаимодействия.
В настоящее время методами численного эксперимента удается
решить широкий круг задач по расчету объемных,
поверхностных и кинетических свойств систем, состоящих из
частиц с молекулярным строением различной сложности: от
сферически симметричных молекул до длинных молекулярных
цепочек и даже полимеров. Трудоемкость методов численного
эксперимента остается достаточно высокой, поэтому пока они
не нашли широкое применение в инженерных методах расчета и
проектирования химико-технологических процессов. Однако
необходимо понимать, что эти методы оказываются
незаменимыми при моделировании наноразмерных объектов. В
настоящее время не существует других теоретических подходов
для разработки нанотехнологий, кроме методов квантовой
механики и молекулярно-динамического моделирования,
поэтому развитие методов численного эксперимента является
весьма перспективным.
Несложно заметить, что основными исходными данными
для проведения расчетов методами МК и МД являются
потенциалы межмолекулярного взаимодействия. Теоретически
рассчитать эти потенциалы можно только методами квантовой
химии. Учитывая ее современное развитие, на практике широко
используют различные модели потенциалов межмолекулярного
взаимодействия.
Модели потенциалов межмолекулярного
взаимодействия
Для адекватного описания поведения системы требуются
значения для потенциалов межмолекулярного взаимодействия
(типичный вид потенциала межмолекулярного взаимодействия
представлен на рис. 4.6).
138
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Теоретически
их
можно
получить
квантовомеханическими расчетами, но ввиду малых значений энергии
взаимодействия молекул, по сравнению с химическими
взаимодействиями, вычисленные этим способом данные лежат в
пределах погрешности расчета, поэтому используют ряд
моделей.
Рис. 4.6. Типичный вид потенциала межмолекулярного
взаимодействия
Модель твердых сфер
Самая простая из всех возможных. Учитывает только
отталкивательную часть взаимодействия.
ϕ = ∞ при r = σ .
ϕ = 0 при r > σ .
Модель степенных потенциалов
Более прогрессивная модель, учитывающая вклад во
взаимодействие как отталкивания, так и притяжения.
ϕ = Ar − n − Br − m .
Частным
случаем
модели
является
потенциал
Леннарда−Джонса:
  σ 12  σ 6 
ϕ = 4ε    −    ,
 r 
 r  

139
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
где σ − эффективный диаметр молекулы, ε − величина
«потенциальной ямы», характеризует энергию максимального
сближения молекул.
Шестая степень при ветке притяжения подтверждена
квантово−механическими расчетами.
Модель Букингема
ϕ = Ar −αr − Br −6 .
Ветка отталкивания задается через экспоненциальную
зависимость.
Модель центр−
−центровых взаимодействий
(атом−
−атомная схема)
Применяется для несферических молекул. Молекула
условно делится на несколько атомных центров, которые можно
считать сферическими. Для каждого из этих центров по одной
из рассмотренных выше моделей находится значение
потенциала межмолекулярного взаимодействия с другими
симметричными молекулами или центрами, выделенными
подобным
же
способом.
Значение
потенциала
межмолекулярного взаимодействия для молекулы в целом
находится суммированием полученных промежуточных
значений.
ϕ12 = ϕα ϕ β .
∑
α ,β
Атом−атомная схема взаимодействий двумя молекулами
СО представлена на рис. 4.7.
Рис. 4.7. Атом-атомная схема взаимодействий между молекул
СО:
140
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
1, 1’ – атомы С; 2, 2’ – атомы О
В заключение необходимо сказать, что МК и МД
являются не единственными фундаментальными методами
молекулярно-статистической теории. В настоящее время
активно развивается теория интегральных уравнений для
частичных функций распределения. Эта теория берет свое
начало с 40-х г. прошлого столетия с работ Боголюбова, Борна,
Грина, Кирквуда, Ивона. Было показано, что теория Гиббса
имеет серьезное продолжение, основанное на том, что
размерность задачи может быть существенно сокращена.
Оказалось, что для описания макросвойств рассматриваемой
системы достаточно знать поведение функций, зависящих от
координат небольшого числа молекул (двух, трех). Эта теория
прекрасно зарекомендовала себя при описании различных
термодинамических свойств простых систем с центральными
межмолекулярными взаимодействиями, и акцент исследований
смещается сейчас на более сложные варианты молекулярных
систем. Ознакомиться с данной теорией, а также более подробно
с методами численного эксперимента можно в следующих
монографиях [17−21].
141
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Заключение
Математическое моделирование в настоящее время
превратилось в универсальный метод, широко используемый в
области химической технологии. Он является мощным
средством синтеза, анализа и интенсификации различных ХТП,
позволяет оценить количественно характеристики процесса,
сравнить их между собою и выбрать наилучшие на стадии
проектирования, когда еще не вложены средства в разработку
промышленного производства и сравнительно легко отказаться
от ошибочно принятого решения.
Насколько обширен круг целей и задач химической
технологии настолько разнообразны и математические описания
ХТП. Рассмотренные в настоящем учебном пособии
математические модели не исчерпывают всего многообразия и
степени полноты математического описания ХТП. Здесь
изложены основные подходы построения математического
описания в процессах, в которых осуществляется перенос тепла,
массы и протекает химическая реакция. Особо выделено то, что
для
получения
результатов
помимо
описания
гидродинамической структуры потока важным является
постановка граничных условий, а также надежное описание
физико-химических свойств рабочих агентов, включая условия
фазового равновесия и коэффициентов переноса субстанций.
Примеры построения расчетных алгоритмов рассмотренных
здесь задач, их программных модулей в пакете Matcad и анализ
решения будут представлены во второй части этого учебного
пособия − в лабораторном практикуме по моделированию ХТП.
142
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Библиографический список
1. Хартман, К. Планирование эксперимента в исследовании
технологических процессов / К. Хартман, Э. Лецкий,
В.Шефер. – М.: Мир, 1977. – 552 с.
2. Разинов, А.И. Теоретические основы процессов химической
технологии: учеб. пособие / А.И. Разинов, О.В. Маминов,
Г.С. Дьяконов. – Казань. – Изд-во Казан. гос. технол. ун-та.
2005. – 362 с.
3. Базаров, И.П. Термодинамика. / И.П. Базаров. – М.:
Высш.шк., 1991. – 376 с.
4. Смирнова, Н.А. Молекулярные теории растворов / Н.А.
Смирнова. – Л.: Химия, 1987. – 336 с.
5. Рид, Р. Свойства газов и жидкостей: справочное пособие / Р.
Рид, Дж. Праусниц, Т. Шервуд; пер. с англ. под ред. Б.И.
Соколова. – 3-е изд., перераб. и доп. – Л.: Химия, 1982. – 592
с.
6. Де Донде, Т., Термодинамическая теория сродства (книга
принципов) / Т. де Донде, П.М. Риссельберг. – М.:
Металлургия, 1984. – 136 с.
7. Разинов, А.И. Явления переноса: учеб. пособие / А.И.
Разинов, Г.С. Дьяконов. – Казань: Изд-во Казан. гос. технол.
ун-та. – Казань, 2002. – 136 с.
8. Кафаров, В.В. Математическое моделирование основных
процессов химических производств: учеб. пособие для вузов
/ В.В. Кафаров, М.Б. Глебов. – М.:Высш.шк., 1991. – 400 с.
9. Павлов, К.Ф. Примеры и задачи по курсу процессов и
аппаратов химической технологии: учеб. пособие для вузов /
К.Ф. Павлов, П.Г. Романков, А.А. Носков. – 13-е изд.,
стереотипное. Перепечатка с издания 1987 г. – Л.: Химия,
2006. – 576с.
10. Плановский, А.Н. Процессы и аппараты химической
технологии / А.Н. Плановский, В.Н. Рамм, С.З. Каган. – 2-е
изд. – М.: Госхимиздат, 1962. – 847 с.
11. Кафаров, В.В. Основы массопередачи / В.В. Кафаров. - М.:
Высш. школа, 1979. – 439 с.
12. Рамм, В.М. Абсорбция газов / В.М. Рамм. – М.: Химия, 1976.
– 654 с.
143
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
13. Yuonglove, B.A. Thermophysical properties of fluids. II.
Methane, Ethane, Propane, Isobutane, and Normal Butane / В.А.
Yuonglove, J.F. Ely. // J. Phys. Chem. Ref. Data, 1987. – V.16. №4. – P.577.
14. Уэйлес, C. Фазовые равновесия в химической технологии / С.
Уэйлес. – М.:Мир, 1989. – 304 с.
15. Хилл, T. Статистическая механика / Т. Хилл. – М.: ИЛ, 1960.
– 435 с.
16. Equation of state calculation by fast computing machines / N.A.
Metropolis, [et. all] // Chem. Phys. – 1953. – V.21 –P. 1087–
1092.
17. Фишер, И.З. Статистическая теория жидкостей / И.З. Фишер.
– М.: Гос.изд.физико-математической лит., 1961.
18. Крокстон, К. Физика жидкого состояния / К. Крокстон; пер. с
англ. – М.: Мир, 1978. – 400 с.
19. Балеску, Р. Равновесная и неравновесная статистическая
механика./ Р. Балеску; пер. с англ. – М.: Мир, 1978.
20. Allen, М.Р. Computer simulation of liquids / М.Р. Allen, D.J.
Tildesley – Oxford: Clarendon press, 1987. – 385 c.
21. Frenkel, D. Understanding molecular simulation. From
algorithms to applications. / D. Frenkel, B. Smit. – Academic
press, 1996. – 638 c.
144
Документ
Категория
Без категории
Просмотров
716
Размер файла
868 Кб
Теги
технологическая, процессов, моделирование, химиков, 158, математические
1/--страниц
Пожаловаться на содержимое документа