О новом подходе к оптимальному проектированию ректификационной колонны Текст научной статьи по специальности «Математика»

УДК 66.048.3

Ф. У. Закирова, Н. Н. Зиятдинов, И. И. Емельянов

О НОВОМ ПОДХОДЕ К ОПТИМАЛЬНОМУ ПРОЕКТИРОВАНИЮ

РЕКТИФИКАЦИОННОЙ КОЛОННЫ

Ключевые слова: оптимальное проектирование, ректификационная колонна, дискретно-непрерывное нелинейное программирование, метод ветвей и границ.

Представлено сравнение авторского подхода к решению задачи оптимального проектирования ректификационной колонны с разработанным раннее подходом. Оценка эффективности подходов проводится по времени расчета.

Keywords: optimal design, distillation column, discrete-continuous nonlinear programming, branch and bounds method.

The paper presents a comparative analysis of two approaches to solving the optimal design problem for a distillation column; one of them is newly developed by the authors of the paper. The approaches are compared and their efficiency is assessed based on the time for calculation.

Введение

Проектирование оптимальных ректификационных колонн с использованием потарелочных математических моделей многокомпонентной ректификации является исключительно сложной проблемой [1]. Трудность решения рассматриваемой задачи обусловлена наличием как непрерывных поисковых переменных, так и дискретных поисковых переменных, характеризующих число тарелок в укрепляющей и исчерпывающей частях колонны, а также сложностью математического описания самого процесса ректификации [2]. В связи с этим рассматриваемая задача формулируется как задача смешанного дискретно-непрерывного программирования. Новый подход к решению этой задачи изложен в работе [3]. В настоящей статье исследуется эффективность рассматриваемого подхода в сравнении со раннее разработанным [4].

Теоретический анализ

В ряде работ для облегчения задачи поиска оптимальной конструкции и режима работы колонны используют упрощенные подходы. Один из них основан на расчете минимального флегмового числа и заключается в следующем. Сначала по уравнению Андервуда [5-7] вычисляется минимальное флегмо-вое число. Далее, оптимальное флегмовое число рассчитывается по формуле:

Яор1 =в* Ят;п , (1)

где в - коэффициент избытка флегмы.

Этот коэффициент является эмпирическим, и в различных работах предлагаются разные способы его получения и интервалы его изменения. Как правило, интервал лежит в пределах в =(1,05 - 1,5) [2, 3].

Также к группе упрощенных методов можно отнести метод Львова-Серафимова, по которому определение числа теоретических тарелок сводится к расчету бинарной смеси из определенных пар компонентов [3].

Однако упрощенные методы не отвечают требованиям точности и оптимальности по технико-экономическому критерию проектирования ректификационных колонн.

Этого недостатка лишены потарелочные математические модели многокомпонентной ректификации, так называемые MESH модели ректификационных колонн. Аббревиатура MESH означает: Mass - уравнения покомпонентного и общего материального баланса; Equilibrium - уравнения равновесия, включающие уравнения точки кипения и точки росы; Summation - уравнения нормировки составов; Heat - теплота, энтальпия или уравнения теплового баланса [1]. На рисунках 1 и 2 представлены блок схемы установки ректификации как сложной химико-технологической системы и тарелки ректификационной установки соответственно.

. ф i

1

2

Л/т-7

I I

Л/т

L

И/ ►

Рис. 1 - Информационная блок схема колонны ректификации

к-ая тарелка

v и^л

ц

Рис. 2 - Информационная блок схема к-ой тарелки колонны ректификации

F

k

F

На данных блок схемах и при дальнейшем описании модели ректификационной колонны индекс k (к = 1,...,NT) будет обозначать номер тарелки в колонне (нумерация сверху вниз), индекс / -

номер компонента смеси (/ = 1.....N0), где КТ -

число тарелок в колонне, N0 - число компонентов в смеси. Мольные доли /-го компонента на к-ой тарелке в жидкой и газовой фазе обозначим через хк и ук соответственно. Тарелки в колонне связаны между собой потоками пара Ук и жидкости Ц. Так же имеется возможность подачи на каждую тарелку потока питания р.

Энтальпия каждого потока (Ну - парового, И/ - жидкостного) рассчитана на основе температуры на тарелках Тк. Давление в колонне определено на каждой тарелке Рк.

Для применения строгих методов необходимо сформулировать математическую модель проектируемой системы, критерий оптимальности, ограничения и указать поисковые переменные, в пространстве которых будет проводиться оптимизация [1, 5, 7].

МЕ8И - математическая модель ректификационной колонны представленной на рисунке 1, включает следующие уравнения:

Уравнение покомпонентного материального

баланса:

Р^к + Цк-1х/,к-1 + Ч+1У/,к-1 - - КУ к = 0, (2) / = 1.....N0, к=1,..,ЛГ.

Уравнения фазового равновесия:

У к = КкХк, (3)

Кк = К(Тк, Рк, Хк ). (4)

Уравнение, связывающее равновесную и рабочую концентрации /-го компонента через Е,к -локальный эффективный коэффициент полезного действия тарелки, имеет вид:

у к = у,,к+1 + Е к (у'к- у ,к+1). (5)

Стехиометрические соотношения для составов жидкой и паровой фаз:

X Хк = 1, X Ук = 1, к = 1,...,ЛГ. (6)

Уравнения теплового баланса: РНРк + ¿.к-Дк-1 + К+Як+1 - 1-кНк - УЛк = о ,(7)

к = 1,...,ЛГ, Нк = Н(Тк ,Рк, Хк), Нк = Н(Тк ,Рк, Ук), Нк = Н(Т, Р ).

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

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

С использованием приведенной математической модели (2)-(7), для проведения расчетов в «поверочной» постановке, в настоящее время разработан ряд эффективных методов, подробный обзор которых дан в [5].

Сложность решения задачи расчета ректификационных колонн в «проектной» постановке, с использованием модели (2)-(7), заключается в том, что наряду с непрерывными переменными, характеризующими режим работы ректификационной колонны, необходимо определить числа тарелок в укрепляющей и исчерпывающей частях колонны, которые являются дискретными переменными.. Таким образом, задача оптимального проектирования ректификационной колонны формулируется как задача смешанного дискретно-непрерывного программирования [1, 8].

В работах [4, 9] приведен подход к оптимальному проектированию систем простых ректификационных колонн, основанный на идеях метода ветвей и границ. Приведем формализованную постановку задачи оптимального проектирования для одной колонны:

* = т ¡П* (х,и,тв), (8)

х ,и ,т$

р (х ,и ,тз) = 0, в = 1,2,

1 < тв < ттах, в в '

у (х ,и ) < 0,

(9) (10)

где

рк (х ,и ,т) =

,ио), в = 1

рДх ,Р ), в = 2

р (х), 1< к< тв, при в=1, (11) вк

1< к < тв, при в = 2

Рт3+1(х 'ит3+1)в = 2

где х, и - переменные состояния и управляющие переменные в ректификационной колонне; в - номер укрепляющей (в = 1) и исчерпывающей (в = 2) секции колонны; тз - число тарелок в укрепляющей (исчерпывающей) секции колонны, которое в общем случае может принимать любые целые значения в пределах от 1 до т™* - заданное максимальное число тарелок в в-ой секции ректификационной колонны; * (х ,и ,тв) - суммарные приведенные капитальные и эксплуатационные затраты колонны; уравнения (9), (11) - математические модели укрепляющей и исчерпывающей секций колонны, р10(х ,и0) - математическая модель дефлегматора в колонне; фзк(х ) - математическая модель к -ой тарелки в колонне; р21(х ,Р ) - математическая модель тарелки питания в колонне;

р (х , и ) - математическая модель кипятиль-тт+1 тт+1

ника в колонне; неравенства (10) - проектные ограничения.

Для простоты дальнейшего изложения в задаче (8) будут опущены уравнения (11).

Для решения задачи (8) в работах [4, 9] предлагается подход к определению нижних оценок критерия оптимальности, позволяющий перейти от дискретных поисковых переменных - числа тарелок в колоннах, к непрерывным структурным параметрам. В основе подхода лежит модификация уравнения, связывающего рабочие и равновесные концентрации компонента через локальный эффективный коэффициент полезного действия тарелки, путем введения структурного параметра. Для этого в каждом уравнении, связывающем рабочие и равновесные концентрации разделяемых компонентов (5)

вводится искусственная переменная азк следующим

образом

VI = УТ +ап1 № - /Г),

(12)

где укз - состав паровой фазы, покидающий к -ую *к

тарелку колонны, уз состав равновесной паровой

фазы, покидающий к -ую тарелку колонны, г/з -локальный эффективный коэффициент полезного действия тарелки. Легко проверить, что если

азк = 1, то к -ая тарелка присутствует в в -ой секции колонны, если азк = 0, то этой тарелки нет. С другой стороны параметры азк могут принимать и непрерывные значения в интервале 0 < азк < 1. Тогда нижняя оценка критерия оптимальности задачи оптимального проектирования простой ректификационной колонны, получается решением следующей задачи:

1 = т1пг(х,и,а3), (13)

х ,и ,а$

р (х , и ,аз) = 0, в = 1,2 , (14)

¥ (х ,и ) < 0, (15)

0 < а3 < 1, а3 = (а31,а32,...,азтГ)г . (16)

Для получения верхней оценки критерия оптимальности для расчета числа тарелок в колонне предложено использовать округленные значения суммы структурных параметров, полученных в результате решения задачи получения нижней оценки.

Проблема, с которой мы здесь сталкиваемся - это большая размерность задачи оптимизации, которая равна сумме структурных переменных. Так например, для решения задачи проектирования одной ректификационной колонны, состоящей из 60 тарелок, число поисковых структурных переменных для поиска нижней оценки нулевой вершины будет равно 60. Понятно, что при решении задачи проектирования системы ректификационных колонн это число будет увеличиваться в арифметической прогрессии и строго зависит от начального приближения числа тарелок в каждой колонне.

Описание предложенного алгоритма

Для сокращения числа поисковых переменных азк был предложен новый подход [10-11]. Для этого на первой итерации в каждой секции колонны

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

Следовательно, на р-ой итерации в каждой секции колонны будет г = 2р-1 групп параметров

азк. В каждой из этих групп параметры азк будут принимать одинаковые значения. Таким образом, на первой итерации мы будем иметь:

а = а = ■■■ = азт1 =а, з=1,2 (17)

На р-ой итерации мы будем иметь:

а = а = ■■■=а = 4^

ав,ц+1 = ав,д+2 = ■■■ = ав.2а = аР2 ,

з,2а ^з2 '

а»,(г-1)а+1 = ав,(г-1)а+2 = ■■■ = авте = аРг , (18)

где а = [тв / 2р-1] - округленное до ближайшего целого число параметров в каждой группе.

Таким образом, на р-ой итерации мы будем решать следующую задачу:

1(р) = т1п 1 (х ,и ,ав), (19)

х ,и ,а8

р (х , и ,а3) = 0, в = 1,2, (20)

¥ (х ,и ) < 0, (21)

0 <азк < 1, аз = (аз1,аз2,...,азтГ )т , (22)

аз1 =аз2 = ■■■ = аза =а3р^

3,2а = аз2 '

з1 з2

з,а+1 = аз,а+2

аз.а+1 = аз.а+2 = ■■■ = аз.2а = ар ,

аз,(г-1)а+1 = аз,(г-1)а+2 = ■■■ = азШе = аРг. (23) В качестве начального приближения, для каждой последующей итерации, используются значения поисковых переменных, полученные на предыдущей итерации. Пусть хр ,ир ,аресть решение

задачи (19). Тогда начальные значения хр0,ир0,ар0 поисковых переменных на р -ой итерации будут связаны с оптимальными значениями, полученными на (р - 1) итерации следующими соотношениями:

хр0 = х( р-1)* ир0 = и(р-1)*,

ар0 = аг.

(24)

(25)

(26)

Рассмотрим решение задачи (19) на первой итерации. Задача (19) отличается от задачи (13) тем, что в ней добавлены ограничения (23). Отсюда оптимальное значение целевой функции, полученное на первой итерации будет заведомо больше оптимального значения целевой функции задачи (13):

1' < 1(1). (27)

На всех последующих итерациях число ограничений (23) будет уменьшаться и, в конце концов,

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

*(р) < *(Р-1). (28) Итерационная процедура заканчивается при выполнении условия:

f (p-1) _ f (p) < £,

(29)

где £ - достаточно малая величина.

Верхняя оценка оптимального значения получается решением следующей задачи:

п1 = min f (x u m), (30)

x ,u

p (x ,u ,as) = 0, s = 1,2, (31)

¥ (x ,u ) < 0, (32)

m™ m™x

где ms = [ £ (askП +1; [ £ (ask)*] - есть округ-k=1 k=1 ление до ближайшего меньшего целого величины

m™

£ (ask )*; (ask) - получены решением задачи k=1

(19).

Постановка и решение задачи оптимального проектирование ректификационной колонны

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

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

M ££ask+Mretr + Mod

F= Pmm*—^---+ FWCL + PvQeb , (33)

где ask - структурный параметр; pmet - цена металла, руб/кг; t - срок окупаемости проектируемой колонны, ч; Meeb - масса металла (кг), расходуемого на изготовление кипятильника колонны; Mcond -масса металла (кг), расходуемого на изготовление

конденсатора колонны; М - масса металла (кг), расходуемого на изготовление одной тарелки и секции кожуха колонны высотой, равной межтарельчатому

расстоянию, рт - цена хладагента, руб/м3; Ош -расход хладагента в колонне, м3/ч; р>у - цена теплоносителя, руб/Гкал; ОгеЬ - тепловая нагрузка на куб колонны, Гкал/ч. М , МееЬ, Моопб - являются функциями конструктивных и режимных параметров колонны.

Расчет эффективности тарелки массообмен-ных устройств является сложной математической задачей и в настоящей статье не рассматривается. В данной работе эффективность тарелок была рассчитана с помощью выражения зависимости от состава питания и его вязкости [12] ц = 0,62. Начальное

приближение числа тарелок для колонны было принято равным 100. Были наложены ограничения на

температуру верхнего потока /-ой колонны, Т^ > 40 °С, при которой возможна конденсация продуктов промышленной водой. Требования на качество следующие: для дистиллята колонны содержание изо-пентана не более 0,0008 массовых долей, для куба колонны содержание н-бутана не более 0,0008 массовых долей.

В качестве поисковых переменных были приняты расходы флегмы в колонне, температура куба колонн, давление верха колонны, переменные

ак.

Таким образом, общая размерность задачи вычисления нижней оценки с помощью уравнений (13)-(16) была равна 103 (первый подход). Размерность задачи вычисления нижней оценки с помощью уравнений (19)-(23) на первой итерации была равна 5, на второй итерации была равна 7 и в общем случае, на ^-ой итерации будет равна 2-1+3 (второй подход). Удвоение числа структурных параметров происходит до тех пор, пока значение нижней оценки на двух последующих итерациях не совпадут с заданной точностью.

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

Ход решения задачи с помощью первого и второго подходов представлены в таблицах 1 и 2 соответственно. Пустые значения, соответствующие верхней оценке какой-либо из вершин, означают, что в этой вершине значение нижней оценки больше наименьшего значения верхней оценки всех существующих вершин и такая вершина исключается из дальнейшего рассмотрения.

Следует отметить, что условие (29) при использовании предложенного подхода выполнялось при р < 3. Из представленных в таблицах результатов видно, что общее время расчета при использовании второго подхода сократилось более чем в 5 раз.

Материальный баланс и характеристики колонны для оптимальной вершины 5.2, полученные с помощью второго подхода, приведены в таблицах 3-4. Следует отметить, что данные результаты сов-

падают с решением, полученным с использованием первого подхода.

Таблица 1 - Решение задачи с использованием

первого подхода (13)

№ вер шины Нижняя оценка критерия Верхняя оценка Время расчета, вр.ед. Точность

mj ml Критерий

0.1 179546 24 24 180932 2.07 0,76

1.1 177166 29 26 178046 3.54 0,49

1.2 180926 2.45

2.1 177091 31 29 177764 1.4 0,38

2.2 177350 29 24 178690 2.78 0,75

3.1 177126 32 28 178036 0.98 0,51

3.2 176976 28 28 177671 1.7 0,39

4.1 178372 29 32 178898 0.62 0,29

4.2 176751 29 26 178046 1 0,72

5.1 177199 28 27 177946 0.52 0,41

5.2 177374 30 28 177446 1.09 0,04

Общее время расчета 18.15

Таблица 2 - Решение задачи с использованием предлагаемого подхода (19)

№ вер ши ны Нижняя оценка критерия Верхняя оценка Время расчета, вр.ед Точность

m m2 Критерий

0.1 177646 29 27 178046 0,22 0,22

1.1 177026 29 27 178046 0,14 0,57

1.2 181826 0,4

2.1 176876 28 27 177952 0,1 0,6

2.2 177766 30 24 178973 0,59 0,67

3.1 177435 32 27 178004 0,18 0,32

3.2 177348 28 26 177942 0,51 0,33

4.1 179362 27 32 180050 0,44 0,38

4.2 177304 28 28 177942 0,47 0,35

5.1 177112 28 27 177942 0,18 0,46

5.2 177426 30 28 177448 0,36 0,01

Общее время расчета 3,58

где т$ - число тарелок в укрепляющей или исчерпывающей секциях ректификационной колонны, нижний индекс - номер секции.

Таблица 3 - Материальный баланс колонны

№ потока/ Параметры потока Сырье Дистиллят Кубовый продукт

Температура, 0С 53,8 40 104

Давление, атм 9 4 4,8

Расход, кг/ч 75010 18550 56460

Состав, масс. доли

изобутанД-С4Н10 0,0464 0,1877 0

бутан, п-С4Н10 0,2012 0,8114 0,0008

изопентан, ьС5Н12 0,1881 0,0008 0,2496

пентан^-СДг 0,1881 0,0001 0,2498

гексан, С6Н14 0,1881 0 0,2499

гептан, С7Н16 0,1881 0 0,2499

Таблица 4 - Оптимальные характеристики колонны

Укрепляющая секция Исчерпывающая секция

Число тарелок 30 Число тарелок 28

Флегмовое число 1,5 Паровое число 1,2

Диаметр секции, м 2 Диаметр секции, м 2

Высота секции, м 18.5 Высота секции, м 17.1

Нагрузка на конденсатор, кДж/ч 1,7*10 7 Нагрузка на ребойлер, кДж/ч 2,3*1 07

Литература

1. L.T. Biegler, I.E. Grossmann, A.W. Westerberg, Systematic methods of chemical process design. New Jersey, 1997, 660 p.

2. Y. D. Land, L.T. Biegler, AIChE Journal, 48, 3, 582-593 (2002).

3. G. M. Ostrovskii, N. N. Ziyatdinov, F. U. Mustafina, Doklady Chemistr, 450, 144-147 (2013).

4. G. M. Ostrovskii, N. N. Ziyatdinov, T. V. Lapteva, and N. Kh. Bogula, Theor. found. chem. eng., 45, 1, 89-98 (2011).

5. Ю.А. Комиссаров, Л.С. Гордеев, Д.П. Вент, Научные основы процессов ректификации: Т.2. Учебное пособие для вузов. Химия, Москва, 2004. 416 с.

6. Александров И.А. Ректификационные и абсорбционные аппараты. Методы расчета и основы конструирования. Химия, Москва, 1978.

7. Н.В. Лисицын, В.К. Викторов, Н.В. Кузичкин, Химико-технологические системы: Оптимизация и ресурсосбережение. Менделеев, Санкт-Петербург, 2007. 312 с.

8. Г.М. Островский, Н.Н.Зиятдинов, Т.В. Лаптева, Д.А. Рыжов, Теор. осн. хим. технологии, 42, 4, 401-412 (2008).

9. Н.Н. Зиятдинов. Н.Ю. Богула, Т.В. Лаптева, Г.М. Островский, Вестник Казанского технологического университета, 5, 118-123 (2011).

10. G. M. Ostrovskii, N. N. Ziyatdinov, F. U. Mustafina, D. A. Rygov, Theor. found. chem. eng., 47, 6, 646-655 (2013)

11. Н.Н. Зиятдинов, Ф.У. Закирова, Д.А.Рыжов, А.А. Караванов, Вестник Казанского технологического университета, 24, 110-118 (2013).

12. O'connell, Trans A.I.Ch.E, 42, 4, 741-755 (1946).

© Ф. У. Закирова - асс. каф. системотехники КНИТУ, firyuza.zakirova@mail.ru; Н. Н. Зиятдинов - д-р техн. наук, проф., зав. каф. системотехники КНИТУ, nnziat@yandex.ru; И. И. Емельянов - асс., асп. той же кафедры, friend41 @mail.ru.