Dynamic model of a waste heat boiler

Full Text

Введение

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

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

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

В работе рассматривается задача реализации динамической модели в пакете MATLAB^© Simulink^©, описывающей поведение котла-утилизатора водогрейного КУВ-24-114 производства ЗАО «Энергомаш (Белгород)-БЗЭМ». Котел предназначен для нагрева сетевой воды за счет использования тепла горячих выхлопных газов газотурбинной установки ГТЭС-16ПА. Котел представляет собой теплообменный аппарат с перекрестным движением потоков воды и выхлопных газов. Система автоматического управления обеспечивает регистрацию температуры и давления выхлопных газов на входе и выходе котла, температуры и расхода воды на входе и температуры на выходе котла. Модель должна обеспечивать расчет температуры воды на выходе при известных остальных параметрах функционирования котла.

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

Часто для моделирования теплообменных аппаратов различных конструкций используется метод конечных объемов в комбинации с другими подходами, обеспечивающими повышение точности. В [6] исследуется точность и скорость моделирования пластинчатого теплообменника в зависимости от количества элементов элементарного объема. В [7] сравнивается поведение моделей теплообменных аппаратов, полученных с помощью метода конечных объемов и с помощью метода движущейся границы. Авторы делают выводы, что метод движущейся границы требует меньше вычислительных ресурсов по сравнению с методом конечных объемов. В [8] проводится исследование термодинамических свойств системы аккумулирования тепла. Для моделирования применяется метод конечных объемов. Полученная модель используется для оптимизации параметров системы.

Для решения задач синтеза автоматических систем управления особенный интерес представляют модели, полученные аналитическими методами [3, 9], т. к. анализ аналитического решения позволяет предложить структуру регулятора. Идентификация таких моделей требует собственного математического аппарата. Так, в работах [10, 11] представлены специальные методы параметрической идентификации объектов, описанных уравнением теплопроводности.

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

Передаточные функции для прямоточных и противоточных теплообменных аппаратов как объектов с распределкнными параметрами, полученные с помощью методов структурной теории распределкнных систем, представлены в качестве примера в монографии А.Г. Бутковского [12]. Представление противоточных теплообменников в виде двусвязной системы сосредоточенных передаточных функций по каналам «температура потоков на входах» – «температура потоков на выходах» получено в работе [13]. В [14] предложено решение, позволяющее реализовать двухсвязную систему в пакете моделирования динамических систем MATLAB^© Simulink^©. Использование аналитических решений, полученных в [12–14], ограничено постоянной скоростью движения потоков. Поэтому отдельный интерес представляет численно-аналитический подход, обеспечивающий разработку адекватных моделей теплообмена потоков энергоносителей, учитывающий параметрические возмущения, такие как изменение скорости потоков, коэффициентов теплопередачи.

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

Более универсальным кажется подход, предложенный В.А. Ковалем в [16], который основан на спектральном представлении распределенных систем: при разложении дифференциального уравнения в частных производных по базису ортонормированных функций уравнение представляется в виде бесконечномерной системы обыкновенных дифференциальных уравнений. Применение спектрального представления при моделировании потока с переменной скоростью исследовалось в [17]. В [18] на базе спектрального представления распределенных систем представлено законченное решение для моделирования процесса нагрева потока с переменной скоростью и его реализация в пакете MATLAB^© Simulink^©. Решение оформлено в виде S-функции [19] и учитывает изменение основных параметров процесса нагрева потока: скорости, коэффициента теплообмена, температуры на входе и температуры нагрева.

Разработанная в [18] S-функция используется ниже для построения Simulink-модели котла-утилизатора.

Математическая модель котла-утилизатора

Поведение температурного распределения потока воды Ɵ(l,t) в трубке теплообменника описывается уравнением

$\frac{\partial \theta \left(l,t\right)}{\partial t}+v\left(t\right)\frac{\partial \theta \left(l,t\right)}{\partial l}+\beta \left(t\right)\theta \left(l,t\right)=\beta \left(t\right)T_{ст}\left(l,t\right),0\le l\le L,t>0$                       (1)

с начальными и граничными условиями

$\theta \left(l\mathrm{,0}\right)={\theta }_0\left(l\right) , \theta \left(\mathrm{0,}t\right)=g\left(t\right)$.                                                                       (2)

Здесь T_cm(l,t) – температура стенки; β(t) – приведенный коэффициент теплообмена, учитывающий величину коэффициента теплопередачи между материалом стенки и водой, площадь конвективного теплообмена и конструктивные особенности теплообменника; v(t) – скорость потока воды; L – общая длина трубки теплообменника; Ɵ₀(l,t) – начальное распределение температуры; g(t) – функция изменения температуры воды на входе теплообменника.

Температура стенки принята одинаковой по всей длине теплообменника. Она определяется исходя из уравнения теплового баланса:

$\frac{d{T}_{ст}\left(t\right)}{dt}=\frac{\left(Q_{г}\left(t\right)-Q_{в}\left(t\right)\right)}{m_{ст}{c}_{ст}}$,                                                                        (3)

$Q_{г}\left(t\right)=c_{г}{F}_{г}\left(t\right)\Delta T_{г}\left(t\right) , \Delta T_{г}\left(t\right)={Т}_{г2}\left(t\right)-{Т}_{г1}\left(t\right)$,                                               (4)

$Q_{в}\left(t\right)=c_{в}{F}_{в}\left(t\right)\left(\theta (L,t)-\theta \left(\mathrm{0,}t\right)\right)$,                                                                   (5)

где Q_г(t) – тепловой поток от газа к стенке теплообменника; Q_в(t) – тепловой поток от стенки к воде; m_cm – масса стенки; c_cm – теплоемкость материала стенки; c_г, c_в – теплоемкости выхлопных газов и воды соответственно; F_г(t), F_в(t)– массовые расходы выхлопных газов и воды; T_г1(t), T_г2(t) – температура выхлопных газов до и после теплообменника.

В [18] для спектрального представления распределения температуры Ɵ₀(l,t) в виде вектора ${\mathbf{Ф}}_{\mathrm{\theta }}={\left[{\mathrm{\phi }}_{\mathrm{\theta }}\left(\mathrm{h},\mathrm{t}\right)\right]}_{\infty \times 1}^{}$ используется базис ортонормированных функций ${\mathbf{Р}}_0\left(l\right)={\left[P\left(h,l\right)\right]}_{\infty \times 1}^{}$:

$\theta (l,t)=$$\sum _{h=1}^{\infty }P(h,l){\phi }_{\theta }(h,t)={\mathbf{\text{Р}}}_0^{\mathrm{Т}}\text{(z)}{\mathbf{\text{Ф}}}_{\text{Θ}}$;                                                         (6)

$P(h,l)=$$\sqrt{2}\sin \left(\frac{(2h-1)\pi l}{2L}\right)$, h = 1,2,.....;                                                 (7)

${\phi }_{\theta }(h,t)$$=\underset{0}{\overset{L}{\int }}\theta (\xi ,t)P(h,\xi )d\xi$, h = 1,2,.....                                                    (8)

Частная производная по пространственной координате δƟ₀(l,t)/δl имеет спектральное представление [16]

$\frac{\partial \theta \left(l,t\right)}{\partial l}={\mathbf{Р}}_1{\mathbf{Ф}}_{\mathrm{\theta }}+{\mathbf{Г}}_0$,                                                                                (9)

где ${\mathbf{Р}}_1={\left[P_1(h,\overline{h})\right]}_{\infty \times \infty }^{}$ – матрица операционного дифференциала первого порядка [16], компоненты которой определяются как

$P_1(h,\overline{h})=$$\underset{0}{\overset{L}{\int }}P(\overline{h},\xi )\left(\frac{\partial P(h,\xi )}{\partial \xi }\right)d\xi$, $h$ = 1,2,....., $\overline{h}$ = 1,2,.....,                               (10)

а вектор ${\mathbf{Г}}_0={\left[{\phi }_g\left(h,t\right)\right]}_{\infty \times 1}^{}$ учитывает ненулевые граничные условия [16].

Уравнение (1) при нулевых граничных условиях g(t) = 0 в спектральном представлении имеет вид

$\frac{\partial {\mathbf{Ф}}_{\theta }}{\partial t}={\mathbf{К}}_1{\mathbf{Р}}_{\mathbf{1}}{\mathbf{Ф}}_{\mathrm{\theta }}+{\mathbf{К}}_2{\mathbf{Ф}}_{\mathrm{\theta }}+{\mathbf{К}}_3\cdot {\mathbf{Ф}}_{\mathit{T}}$,                                                                     (11)

где К₁ =-v(t)Е, К₂ =-β(t)Е, К₃ =-β(t)Е, E – единичная диагональная матрица. Компоненты матрицы Ф_Т = ${\left[{\phi }_T(h,t)\right]}_{\infty \times 1}^{}$ вычисляются для T_cm(l,t) аналогично (8).

Система линейных уравнений (11) используется для реализации модели пространства состояний температуры на выходе теплообменника. Матрицы представления в пространстве состояний рассчитываются на базе ограниченного первыми N членами бесконечного ряда (7).

Учет ненулевых граничных условий выполнен за счет участка дополнительного нагрева в начальной части нагревателя [0,d] с температурой нагрева T_g. Выражение для расчета температуры T_g получено из уравнения для установившегося режима следующим образом. Если рассмотреть режим при значениях параметров β(t) = β_c = const, v(t) = v_c = const, g(t) = g_c = const и нулевой температуре стенки, T_cm(l,t) = 0, то

$\theta \left(l,\infty \right)=g_c\exp \left(-{\beta }_{c}l/v_c\right)$.                                                                               (12)

Для случая, когда температура на входе равна нулю g(t) = 0, а температура стенки задана выражением

$T_{ст}\left(l,t\right)$$=\left\{\begin{array}{l}{T}_g,\text{\hspace{1em}при\hspace{1em}}l\le d;\\ \mathrm{0,}\text{\hspace{1em}при\hspace{1em}}l>d,\end{array}\right.$                                                                           (13)

температура на участке l > d в установившемся режиме будет определяться выражением

$\theta (l,\infty )=T_g\exp \left(-{\beta }_{c}l/v_c\right)\left(\exp \left(-{\beta }_{c}d/v_c\right)-1\right)$, l > d.                                    (14)

Из равенства (12) и (14) получено выражение для расчета температуры на участке дополнительного нагрева:

$T_g=g_c\frac{1}{\exp \left(-{\beta }_{c}d/v_c\right)-1}$.                                                                        (15)

При реализации S-функции принято допущение, что температура на участке дополнительного нагрева для переходных режимов рассчитывается также согласно выражению (15):

$T_g\left(t\right)=g\left(t\right)$$\frac{1}{\exp \left(-\beta \left(t\right)d/v\left(t\right)\right)-1}$.                                                                (16)

Реализация и исследование модели в Matlab® Simulink®

Для применения S-функции sfFlowHeating из [18] при создании модели котла-утилизатора необходимо реализовать расчет температуры стенки, а также приведенной скорости потока, поскольку S-функция sfFlowHeating реализована для теплообменника с приведенной длиной трубки L_пр = 1.

Масса стенки теплообменника m_cm может быть рассчитана исходя из характеристик теплообменника котла-утилизатора (см. таблицу).

 

Характеристики теплообменника котла-утилизатора

Характеристика	Значение
Трубы:
– диаметр, мм – толщина стенки, мм – материал	38 3 сталь марки 20К-18
Расположение труб	шахматный порядок
Поперечный шаг труб, мм	110
Продольный шаг труб, мм	75
Количество труб в одном ряду	20
Число рядов в одном ходе	4
Количество ходов	6
Длина хода, мм	5000
Оребрение труб:
– высота, мм – толщина, мм – шаг, мм	16 1 4

 

По данным таблицы рассчитана масса стали теплообменника m_cm = 19 500 кг, внутренний объем теплообменника (по воде) – 2,123 м³. Площадь поверхности теплообмена по воде составляет 241,6 м². Площадь поверхности теплообмена по выхлопным газам – 3544 м².

Теплоемкость уходящих газов c_г принята равной 1,074 кДж/(кг·ºС) [20].

Теплоемкость стали c_см – 0,46 кДж/(кг·ºС).

Теплоемкость воды c_в – 4,2 кДж/(кг·ºС).

Разработанная Simulink-модель котла-утилизатора представлена рис. 1. Температура стенки рассчитывается с помощью блока-интегратора, на вход которого подается разность тепловых потоков, согласно (3). Поток тепла от выхлопных газов рассчитывается при постоянном массовом расходе F_г(t) = 182 т/час, при изменяющейся разности температур ∆T_г(t). Поток тепла в воду рассчитывается исходя из данных по расходу воды на входе котла-утилизатора, разности температур на входе и на выходе теплообменника. На вход S-функции sf FlowHeating подаются (сверху вниз на рис. 1): температура воды на входе теплообменника, g(t); температура стенки теплообменника, T_cm(t); приведенная скорость потока воды, v_пр(t); приведенный коэффициент теплообмена β(t). Приведенная скорость потока воды рассчитывается согласно выражению

$v_{пр}\left(t\right)=$$\frac{v\left(t\right)}{L}=\frac{1}{L}\cdot \frac{{\rho }_{в}{F}_{в}\left(t\right)}{S}=\frac{{\rho }_{в}}{V_{ТО}}$$F_{в}\left(t\right)$,                                                          (17)

где S – суммарная площадь сечения всех трубок теплообменника; ρ_в– плотность воды; V_ТО – внутренний объем теплообменника. Физический смысл приведенной скорости потока – это величина, обратная транспортному запаздыванию.

 

Рис. 1. Simulink-модель котла-утилизатора

 

Для определения значения приведенного коэффициента теплообмена = const (переменная betta) использовались реальные данные, которые подавались на входы Simulink-модели: перепад температуры выхлопных газов ∆T_г(t) (переменная dTug), массовый расход воды F_в(t) (переменная F_в), температура воды на входе теплообменника g(t) (переменная Tw₁). Значение коэффициента получено в результате решения задачи минимизации:

$v_{пр}\left(t\right)=\frac{v\left(t\right)}{L}=\frac{1}{L}\cdot \frac{{\rho }_{в}{F}_{в}\left(t\right)}{S}=\frac{{\rho }_{в}}{V_{ТО}}{F}_{в}\left(t\right)$,                                                    (18)

Где T_в2(t) – реальные данные о температуре воды на выходе котла-утилизатора (переменная Tw₂).

Приведенный коэффициент теплообмена равен β(t) = 0,00528 c^–1, что соответствует коэффициенту конвективного теплообмена «стенка – вода» α_в = 354,8 Вт/(м²·ºС). Коэффициент конвективного теплообмена между выхлопными газами и поверхностью нагрева исходя из величины теплового потока и площади поверхности теплообмена равен α_в = 30.6 Вт/(м²·ºС). Оба значения согласуются с оценками [21].

 

Рис. 2. Изменение технологических параметров в процессе эксплуатации котла-утилизатора

 

На рис. 2 приведены данные, полученные в процессе эксплуатации котла-утилизатора. Данные использовались для идентификации модели. Газотурбинная установка, от которой поступают выхлопные газы на вход котла-утилизатора, работает в номинальном режиме с максимальной нагрузкой. Расход выхлопных газов постоянен. Температура воды на входе практически не изменяется, расход воды изменяется в существенном диапазоне от 400 до 720 т/час, что соответствует изменению приведенной скорости потока от 0,053 до 0,093 1/с или изменению транспортного запаздывания в диапазоне от 10,8 до 18,9 с. Анализ графиков позволяет заметить, что температура воды на выходе существенно зависит от расхода воды.

 

Рис. 3. Поведение модели котла-утилизатора

 

На рис. 3 приведены результаты моделирования работы котла-утилизатора при полученном в процессе идентификации значении коэффициента теплообмена. При исследовании Simulink-модели (см. рис. 1) настроечные параметры S‑функции sfFlowHeating имели следующие значения: N = 10, d = 0,4. Модель обладает высокой точностью: средняя абсолютная ошибка (MAE) составляет 0,2 °С, что соизмеримо с классами допуска современных термометров сопротивления. Максимальное отклонение модели от реальных данных составляет 1,9 °С.

Заключение

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

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

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

About the authors

Ivan A. Danilushkin

Samara State Technical University

Email: idanilushkin@mail.ru

Ph.D. (Techn.), Associate Professor

Russian Federation, 244, Molodogvardeyskaya str., Samara, 443100

Sergey A. Kolpashchikov

Samara State Technical University

Email: skolpaschikov@mail.ru

Ph.D. (Techn.), Head of Department

Russian Federation, 244, Molodogvardeyskaya str., Samara, 443100

Dmitry S. Lusenko

Samara State Technical University

Author for correspondence.
Email: Dmitry-Lusenko@yandex.ru

Postgraduate Student

Russian Federation, 244, Molodogvardeyskaya str., Samara, 443100

References

Supplementary files