Selection of numerical method for solving ordinary differential equation systems for a high-speed model of hydrocarbons steam cracking

Full Text

Введение

Легкие олефины являются основным сырьем для глубоких переделов нефтехимической и химической промышленности. Продукты пиролиза углеводородов являются сырьем для производства пластиков, синтетических волокон, каучуков и др. [1] и товаров народного потребления и промышленного назначения. В связи с постоянно растущим спросом на легкие олефины для дальнейших переделов присутствует необходимость повышения эффективности пиролизных производств. Одним из способов оптимизации работы печей пиролиза является применение систем усовершенствованного управления технологическим процессом (СУУТП) и оптимизации в реальном времени (RTO – Real-Time Optimization). Для работы таких систем требуется модель процесса пиролиза, не только адекватно предсказывающая основные параметры технологического режима, но и способная также выполнять расчеты оптимизационных сценариев за короткое время, что необходимо для компенсации возникающих возмущений процесса в реальном времени.

Проблематика моделирования и методы исследования

Одной из особенностей процесса пиролиза является одновременное протекание большого количества реакций (до 3288 реакций со 128 компонентами и 20 радикалами [2]). Такое количество реакций и вовлеченных в них компонентов накладывают ограничения на применяемые численные методы решения систем ОДУ, т. к. увеличивается общее количество решаемых уравнений. Для снижения количества расчетных компонентов прибегают к использованию псевдокомпонентов, объединяющих близкие по свойствам вещества [3, 4]. Отдельно стоит вопрос о выборе реакций, которые необходимо включать в решаемую систему для получения требуемой точности расчета, и об определении кинетических параметров используемых реакций.

Решение данной системы осложняется жесткостью системы ОДУ материального баланса, которая связана с большими различиями в кинетических параметрах различных реакций (молекулярных и радикальных), а также тем, что малые концентрации радикалов в потоке (порядка до 10^–30 мольн. долей) вносят существенный вклад в скорости образования/исчезновения молекулярных веществ за счет высоких значений констант скоростей реакций, в которых участвуют радикалы [5].

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

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

Разработка модели пиролиза

Предметом исследования является математическая модель реакционного змеевика печи пиролиза, реализованная в виде программного кода на языке Python.

Разрабатываемая модель состоит из семи пользовательских модулей. Структура модели показана на рис. 1.

В модуле основных классов описаны методы расчета термодинамических свойств индивидуальных компонентов (теплоемкость по уравнению DIPPR 107 [6] и энтальпия интегрированием от энтальпии образования в стандартных условиях) и их смесей. На каждом шаге расчета изменяется состав реакционной смеси и, соответственно, плотность, и объемный расход смеси. Для расчета термодинамических свойств потока на данном этапе используется уравнение состояния идеального газа. Результатом работы модуля основных классов являются рассчитанные термодинамические свойства чистых компонентов и затем многокомпонентного потока на каждом шаге интегрирования.

В модуле кинетики выполняется расчет скорости химических реакций и непосредственно решение реактора идеального вытеснения в стационарной постановке задачи ( $\raisebox{1ex}{$\partial C_i$}\!\left/ \!\raisebox{-1ex}{$\partial t$}\right.=0)$, проверяется материальный баланс реактора. Изменение концентраций реагирующих компонентов описывается следующим уравнением:

$u\frac{\partial C_i}{\partial l}={\sum }^{\text{}}{W}_{ij}$ ,

где $C_i$ – концентрация соответствующего i-го вещества, моль/м³; u – линейная скорость потока, м/с; $\partial l$ – элементарная длина ячейки интегрирования, м; t– длина змеевика, м; $W_i$ – скорость j-й реакции по i-му веществу, моль/(м³∙с).

Скорость реакции определяется по закону действующих масс, константа скорости реакции – по уравнению Аррениуса. Обратные реакции учитываются как самостоятельные реакции.

$W_j=k\cdot \Pi C_i ,$

где $C_i$ – концентрация i-го реагента, моль/м³; k – константа скорости реакции

$k=A\cdot e^{-\frac{E_a}{RT}},$

где A – предэкспоненциальный множитель, с^–1 для одномолекулярных реакций или м³/(моль∙с) для реакций с двумя реагентами; $E_a$ – энергия активации, Дж/моль.

 

Рис. 1. Структура разрабатываемой модели

Fig. 1. Model structure

 

Тепловой баланс на каждом шаге рассчитывается по уравнению:

$\frac{\text{dT}}{d\tau }=\frac{{\displaystyle {\sum }^{\text{}}-\Delta H_j\cdot W_j\cdot R}}{P\cdot C_p},$

где $\Delta H_j$– тепловой эффект j-й реакции, Дж/моль; R – газовая постоянная, (м³∙Мпа)/(моль∙К); P – давление, МПа; $C_p$– теплоемкость реакционной смеси, Дж/(моль∙К)

При этом расчет теплообмена с дымовыми газами печи не реализован, однако доступно задание распределенной по длине змеевика тепловой нагрузки. Введение внешней тепловой нагрузки учитывается в расчете добавлением дополнительного слагаемого:

$\frac{\text{dT}}{d\tau }=\frac{\left({\displaystyle {\sum }^{\text{}}-\Delta H_j\cdot W_j+Q}\right)\cdot R}{P\cdot C_p},$

где Q – удельный тепловой поток к змеевику от дымовых газов, Дж/(с∙м³).

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

При относительно больших шагах интегрирования (от 10^–6 м, в зависимости от условий протекания реакций) расход вещества, рассчитанный с учетом параллельного протекания реакций, в которых в качестве реагентов участвуют одни и те же вещества, может оказаться больше, чем содержание данного вещества в потоке

${\sum }_j{W_i}_j\cdot dv>F_i,$

где dv – объем элементарной ячейки интегрирования, м³; $F_i$ – мольный расход i‑го компонента.

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

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

Представленный набор уравнений численно решается на каждом шаге в модуле интегратора. В качестве граничных условий расчета используется следующий набор ограничений:

$T\left(l_0\right)=T_0$,

$C_i\left(l_0\right)=C_{i0}$,

$P\left(l\right)=\text{const}=P_0$.

В интеграторе реализовано два метода решения ОДУ: явный метод Рунге–Кутты–Фельберга пятого порядка [7] и явный метод Рунге–Кутты–Мерсона пятого порядка с адаптивным шагом [8]. На первых этапах разработки модели использовался только метод со статичным шагом интегрирования. На простых кинетических схемах процесса (до 7–10 реакций) и для коротких змеевиков (до 2 м), на которых проводилась отладка модели, использование неизменного шага интегрирование не приводило к каким-либо проблемам при выполнении расчетов и проведении численных экспериментов. Однако с увеличением количества реакций в используемой кинетической модели и при переходе к параметрам змеевика, близким к промышленным аппаратам, время каждого расчета увеличилось до неприемлемых значений (8 и более часов). Такое поведение расчетной модели, несмотря на приемлемую точность расчета, значительно усложняет отладку модели, а также практически полностью исключает возможность использования разрабатываемой модели для решения задач в промышленности. Внедрение метода решения систем ОДУ с автоматическим подбором величины шага интегрирования позволило решить часть возникших проблем.

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

$k_0=h_0\cdot f\left(x_i,y_i\right),$

$k_1=h_0\cdot f\left(x_i+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.\cdot h_0,y_i+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.\cdot k_0\right),$

$k_1=h_0\cdot f\left(x_i+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.\cdot h_0,y_i+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.\cdot k_0\right),$

$k_3=h_0\cdot f\left(x_i+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.\cdot h_0,y_i+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$8$}\right.\cdot k_0+\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$8$}\right.\cdot k_2\right),$

$k_3=h_0\cdot f\left(x_i+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.\cdot h_0,y_i+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$8$}\right.\cdot k_0+\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$8$}\right.\cdot k_2\right),$

${}_{i+1} =y_i+\frac{k_0+4*k_3+k_4}{6},$

где y_i – значение искомой функции на в точке x_i ; y_i+1– значение функции в точке x_i + h; h – шаг интегрирования.

 

Рис. 2. Используемый алгоритм выбора шага интегрирования

Fig. 2. Step selection algorithm used

 

Для контроля точности на каждом шаге интегрирования вычисляется локальная ошибка округления по формуле:

${}_{i+1} =y_i+\frac{k_0+4*k_3+k_4}{6},$.

Алгоритм выбора шага интегрирования на каждой итерации приведен на рис. 2. Алгоритм основан на методе Рунге–Кутты–Мерсона [8, 9] и упрощен в части расчета величины нового шага.

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

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

Возможность подключения создаваемых пользователем баз данных компонентов и реакций позволяет выполнять расчеты с использованием различных кинетических моделей.

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

Описанная модель была опробована на одном наборе данных действующего производства пиролиза СИБУР-Кстово. Для опробования модели были взяты данные по режиму работы и составам сырья и пирогаза печи пиролиза F-01A за 31.03.2011 г. В расчетах использовалась упрощенная кинетическая модель пиролиза этана, состоящая из 56 реакций по [10]. Данная модель включает в себя основные реакции, свойственные для пиролиза легкого сырья – пропан-бутановой фракции – и применима для расчета печей, работающих на этановых режимах. Параметры некоторых реакций, используемых в модели, приведены в табл. 1. Значительное влияние на работу змеевика оказывает образование отложений кокса на стенках труб [11, 12], однако на данном этапе предлагаемые в [10] реакции коксообразования не учтены в общей схеме превращений для упрощения расчетов.

 

Таблица 1. Используемая кинетическая модель

Table 1. Kinetic model used

Реакция/Reaction	A	Ea	Реакция/Reaction	A	Ea
	[s-1] [ml/mol/s]	[kJ/mol]		[s-1] [ml/mol/s]	[kJ/mol]
C2H6→CH3*+CH3*	4,0∙1016	365,84	C2H4+H*→C2H3*+H2	8,0∙1011	16,71
1-C4H8→C3H5*+CH3*	8,0∙1016	309,3	C2H6+H*→C2H5*+H2	1,0∙1014	40,57
n-C4H10→C2H5*+C2H5*	1,5∙1016	343,39	C2H4+CH3*→C2H3*+CH4	1,0∙1013	54,38
C2H3*+H*→C2H4	1,0∙1013	0	CH3*+CH3*→C2H6	1,3∙1013	0

 

При опробовании модели сравнивались расчетные значения выходов основных продуктов пиролиза этана. Использованные исходные данные и результаты расчетов приведены в табл. 2, 3.

 

Таблица 2. Исходные данные для верификации расчетной модели

Table 2. Initial data for model verification

Параметр Parameter	Значение Value
Расход сырья, кг/ч/Feed rate, kg/hr	8360
Состав сырья, мольн. Доли/Feed composition, mol. fract.
Метан/Methane	0,0313
Этан/Ethane	0,945
Пропан/Propane	0,0019
Пропилен/Propylene	0,0001
Изобутан/i-Butane	0,0001
н-Бутан/n-Butane	0,0037
Бутен/Butene	0,0092
н-Пентан/Pentane	0,0008
Изопентан/i-Pentane	0,0062
Гексан/Hexane	0,0018
Расход пара, кг/ч/Dilution steam rate, kg/hr	3850
Давление, МПа/Pressure, MPa	0,518
Температура/Temperature, K	1124
Длина змеевика, мм/Coil length, mm	7700
Труба змеевика, мм´мм/Coil dimensions, mm×mm	90×6,5
Кол-во змеевиков, шт/Number of coils, pcs.	16

 

Таблица 3. Результаты опробования расчетной модели

Table 3. Model verification results

Состав продукта, мольн. доли Composition, mol. fract.	Факт Actual	Расчет Calculated	Отн. ошибка Rel error, %
Водород/Hydrogene	0,3750	0,3618	–3,5
Метан/Methane	0,0942	0,0573	–39,2
Этан/Ethane	0,1784	0,1775	–0,5
Этилен/Ethylene	0,3295	0,3756	14,0
Ацетилен/Acetylene	0,0041	0,0129	218,3
Пропан/Propane	0,0007	0,0012	67,5
Пропилен/Propylene	0,0065	0,0009	–86,8

 

В текущей версии расчетной модели средняя ошибка по основным компонентам реакционной смеси пиролиза этана: этилену и этану (выделены шрифтом в табл. 3), составляет 7,2 %. По другим компонентам наблюдаются значительные расхождения, связанные, в первую очередь, с несовершенством кинетической модели. Выбор реакций, включаемых в модель, а также кинетические параметры реакций определяют предсказывающую способность модели [4, 13–17]. Доработка используемой кинетической модели является основным фокусом будущей работы.

Экспериментальная часть исследования

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

• Расчет 1 – фиксированный шаг 10^–5 м;
• Расчет 2 – адаптивный шаг.

Расчеты выполнялись на тех же исходных данных, на которых проверялась адекватность модели.

На рис. 3 приведен профиль концентраций для расчета 1.

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

  

Рис. 3. Профили концентраций – Расчет 1

Fig. 3. Concentrations profiles – Case 1

 

Рис. 4. Сравнение результатов расчетов – Профиль концентраций этилена – Расчеты 1 и 2

Fig. 4. Results comparison – Ethylene concentrations profile – Cases 1 and 2

 

Рис. 5. Профили концентраций некоторых радикалов – Расчет 2

Fig. 5. Selected radical species concentrations profiles – Case 2

 

Рис. 6. Профили некоторых скоростей реакций – Расчет 2

Fig. 6. Selected reactions rates profiles – Case 2

 

При этом форма профилей концентраций компонентов соответствует общепринятому представлению о кинетике процесса пиролиза. Это хорошо видно на примере профилей концентраций радикальных компонентов (рис. 5) – явно прослеживаются зоны инициации реакций, продолжения реакций и обрыва цепи реакций.

Однако машинное время, затраченное на выполнение двух расчетов на одной и той же рабочей станции, различается радикально: 11 ч 32 мин 19,1 с для Расчета 1 против 33 мин 49,1 с для Расчета 2. Соответственно число выполненных итераций составляет 770 тыс. для Расчета 1 и 36,2 тыс. для Расчета 2.

На рис. 7 приведены профили локальной ошибки округления для двух расчетов.

В данном случае максимально допустимым значением ошибки округления было принято 10^–11 [18]. Видно, что для некоторых участков расчеты выполняются с точностью округления по нижней границе допустимого диапазона, что связано с необходимостью проверки материального баланса – при допустимом значении оценки локальной ошибки округления шаг интегрирования может быть таким, что материальный баланс в ячейке сохраняться не будет (например, расход компонента в реакции будет превышать его количество в ячейке).

 

Рис. 7. Сравнение результатов расчетов – Ошибка округления – Расчеты 1 и 2

Fig. 7. Results comparison – Local truncation error – Cases 1 and 2

 

Рис. 8. Сравнение результатов расчетов – Шаг интегрирования – Расчеты 1 и 2

Fig. 8. Results comparison – Step size – Cases 1 and 2

 

Такое поведение связано с особенностями метода [8] оценки ошибки округления, который не гарантирует достаточно малый шаг для успешного решения системы ОДУ изменения концентраций компонентов. В таких случаях после проверки материального баланса ячейки, если расчетное уменьшение количества вещества превышает его содержание в потоке, шаг интегрирования дополнительно уменьшается.

Однако в некоторых частях змеевика допустимо выполнять расчеты с максимальной ошибкой округления, что позволяет увеличить шаг интегрирования. Так, для Расчета 2 среднее значение оценки локальной ошибки округления составляет 0,52∙10^–11 и 12 % итераций выполнены с ошибкой округления, превышающей среднее значение. При этом все итерации выполнены с ошибкой округления большей, чем при Расчете 1, в котором алгоритм не имеет возможности управлять точностью вычислений, что говорит об избыточной точности расчета на некоторых участках змеевика.

На рис. 8 приведены профили размера шага интегрирования для двух расчетов. Пульсации величины шага для Расчета 2 связаны с описанным выше несовершенством метода оценки локальной ошибки округления.

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

При этом из рис. 8 также видно, что минимальный шаг, который выбирает алгоритм, на некоторых участках составляет порядка 2∙10^–5 м, что не позволяет сократить время расчета при использовании фиксированного шага простым увеличением шага интегрирования. Так, расчет данной задачи с шагом интегрирования 10^–4 м не позволяет свести материальный баланс расчетной ячейки уже после 0,3 м змеевика.

Выводы

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

 

Таблица 4. Сравнение расчетов

Table 4. Results comparison

Параметр Parameter	Расчет 1 Case 1	Расчет 2 Case 2	Различие Difference, %
Кол-во итераций Number of iterations	770000	36214	–95,3
Время расчета Calculation time	11 ч 32 мин 19,1 с	33 мин 49,1 с	–95,1
Состав продуктов (с учетом воды), мольн. доли Products composition (wet), mol. fract.
Вода/Water	0,3228	0,3228	0,00
Водород/Hydrogene	0,2447	0,2450	0,15
Метан/Methane	0,0391	0,0388	–0,92
Этан/Ethane	0,1203	0,1202	–0,03
Этилен/Ethylene	0,2547	0,2543	–0,13
Ацетилен/Acetylene	0,0088	0,0088	–0,31
Пропан/Propane	0,0008	0,0008	–0,05
Пропилен/Propylene	0,0006	0,0006	0,88

 

В работе показано, что использование численных методов решения систем ОДУ при расчетах реакционного змеевика печи пиролиза позволяет значительно сократить затраты машинного времени на расчет при сохранении удовлетворительной точности. Например, при расчете упрощенной схемы пиролиза этана применение адаптивного шага интегрирования позволяет сократить количество расчетных итераций более чем в 21 раз и затрачиваемое на расчет время – в 20 раз. При этом различие по основным компонентам в результатах двух расчетов составляет в среднем 0,31 %. (табл. 4). Сокращение времени расчета достигается за счет увеличения шага интегрирования на тех участках змеевика, для которых характерны пониженные скорости реакций.

Таким образом, для решения промышленных задач, для которых требуется частое (например, 1 раз в 4 часа) обновление результатов расчетов, целесообразно использовать методы интегрирования с адаптивным шагом.

Дальнейшее ускорение расчетов возможно за счет реализации пятишагового неявного метода Гира для решения системы дифференциальных уравнений, который характеризуется высокой стабильностью и точностью при решении жестких систем ОДУ [19, 20].

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

About the authors

Vladimir V. Kozlov

National Research Tomsk Polytechnic University; LLC «NIOST»

Author for correspondence.
Email: kozlovvv@tpu.ru
ORCID iD: 0009-0000-6979-5153

Postgraduate Student, Lead Engineer

Russian Federation, Tomsk; Tomsk

Igor М. Dolganov

National Research Tomsk Polytechnic University

Email: dolganovim@tpu.ru
ORCID iD: 0000-0002-9848-4933

Cand. Sc., Associate Professor

Russian Federation, Tomsk

Stepan S. Slobodin

LLC «SIBUR-Kstovo»

Email: slobodinss@sk.sibur.ru

Process and Technology Innovations Manager

Russian Federation, Kstovo

References

Supplementary files