Численно-аналитическая модель температуры металла барабана парового котла

Полный текст

Введение

Барабаны паровых котлов высокого давления, с рабочим давлением более 10 МПа, представляют собой горизонтально расположенный цилиндр диаметром до 1800 мм изготовленный из листовой легированной стали толщиной 90 мм и более. Регламент растопки котлов ограничивает перепад температуры верхней образующей и нижней образующей барабана, перепад температур насыщения пара и стенки барабана, а также максимальную скорость роста температуры стенки [1]. Перепад и скорость роста температур в контрольных точках стенок барабана входят в список параметров, определяющих действия машиниста котла при растопке. Поэтому, при разработке компьютерного тренажёра машиниста парового котла, возникла задача получения модели, с достаточной точностью описывающей поведение температурного распределения стенки барабана.

В процессе функционирования парового котла в нижней части барабана располагается пароводяная смесь, в верхней – насыщенный пар. Нагрев стенки барабана осуществляется с внутренней стороны от конвективного теплообмена с паром/пароводяной смесью. С наружной стороны барабан теплоизолирован. В процессе растопки котла характеристики теплоносителя существенно изменяются: температура повышается от 30 °C до 340 °C, давление возрастает от 0.1 до 14.5 МПа. Это приводит к изменению параметров теплообмена с внутренней стороны стенки барабана и должно быть учтено при моделировании процесса.

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

Известные программные пакеты компьютерного моделирования, такие как COMSOL Multiphysics^®, ANSYS, Altair Flux™, реализуют решение задач математической физики методом конечных элементов. Все они обладают высокой вычислительной ресурсоёмкостью, кроме того, зачастую требуются дополнительные решения по интеграции полученной модели с программными продуктами сторонних производителей [5].

Если численные модели разрабатываются для решения задач идентификации процессов технологической теплофизики, синтеза систем автоматического управления, систем имитационного моделирования, то необходимо снизить вычислительную ресурсоёмкость алгоритмов моделирования, а также иметь возможность простой реализации моделей в пакетах моделирования динамических систем, таких как Matlab^® Simulink^®, SimInTech, VisSim. В этом случае авторы чаще всего прибегают к методу конечных разностей. Например, в работах [6, 7], посвящённых моделированию температурного распределения стенки барабана котла, авторы используют сетку по пространственной координате, на базе которой переходят к системе линейных дифференциальных уравнений для температур в узловых точках. Такой же подход применён в [8], при решении подобной задачи с переменным коэффициентом теплопередачи.

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

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

Постановка задачи

Температурное распределение по толщине бесконечной пластины с конвективным теплообменом на поверхностях, описывается одномерным параболическим уравнением с граничными условиями третьего рода [2]:

$\frac{\partial Q(x,t)}{\partial t}=a\left(t\right)\frac{{\partial }^{2}Q(x,t)}{\partial x^2}$, $0<x<L$, $t>0$; (1)

$Q\left(x\mathrm{,0}\right)=Q_0\left(x\right)$; (2)

$\left\{\begin{array}{l}\lambda \left(t\right)\frac{\partial Q\left(\mathrm{0,}t\right)}{\partial x}={\alpha }_1\left(t\right)\left(Q_1\left(t\right)-Q\left(\mathrm{0,}t\right)\right),\\ -\lambda \left(t\right)\frac{\partial Q(L,t)}{\partial x}={\alpha }_2\left(t\right)\left(Q_2\left(t\right)-Q(L,t)\right),\end{array}\right.$ (3)

где $Q(x,t)$ – температурное поле пластины; $Q_1(x,t)$, $Q_2(x,t)$ – температуры сред на границах пластины при $x=0$ и $x=L$; $a\left(t\right)=\lambda \left(t\right)/\left(c\left(t\right)\gamma \left(t\right)\right)$ – коэффициент температуропроводности: $\lambda \left(t\right)$ – коэффициент теплопроводности, $c\left(t\right)$ – удельная теплоемкость, $\gamma \left(t\right)$ – плотность материала пластины; ${\alpha }_1\left(t\right)$, ${\alpha }_2\left(t\right)$ – коэффициенты теплопередачи на границах пластины; L – толщина.

При постоянных значениях всех коэффициентов задачи: ${\alpha }_1\left(t\right)={\alpha }_{1C}$, ${\alpha }_2\left(t\right)={\alpha }_{2C}$, $\lambda \left(t\right)={\lambda }_C$, $c\left(t\right)=c_C$, $\gamma \left(t\right)={\gamma }_C$ решение краевой задачи (1) – (3), может быть найдено в терминах структурной теории распределённых систем [16, 17] с помощью передаточной функции объекта с распределёнными параметрами [18]. Стандартизирующая функция для (1) – (3) имеет вид [18]:

$\omega (x,t)=Q_0\left(x\right)\delta \left(t\right)+a_C{g}_2\left(t\right)\delta (x-L)-a_C{g}_1\left(t\right)\delta \left(x\right)$, (4)

где

$g_1\left(t\right)=b_1\left(Q_1\left(t\right)-Q\left(\mathrm{0,}t\right)\right)$,

$g_2\left(t\right)=b_2\left(Q_2\left(t\right)-Q(L,t)\right)$,

$b_1={\alpha }_{1C}/{\lambda }_C$,

$b_2={\alpha }_{2C}/{\lambda }_C$.

Передаточная функция определяется выражением:

$W(x,\xi ,p)=\sum _k\frac{{\varphi }^{*}({\mu }_k,x){\varphi }^{*}({\mu }_k,\xi )}{{||{\varphi }^{*}({\mu }_k,x)||}^2\left(p+a{\mu }_k^2\right)}$, $k=\mathrm{1,2,...}$, (5)

где p – оператор преобразования Лапласа,

${\varphi }^{*}({\mu }_k,x)=\cos {\mu }_{k}x+b_1\frac{\sin {\mu }_{k}x}{{\mu }_k}$,

${||{\varphi }^{*}\left({\mu }_k\right)||}^2=\frac{b_2}{2{\mu }_k^2}\frac{{\mu }_k^2+b_1^2}{{\mu }_k^2+b_2^2}+\frac{b_1}{2{\mu }_k^2}+\frac{L}{2}\left(1+\frac{b_1^2}{{\mu }_k^2}\right)$,

${\mu }_k$ – положительные корни трансцендентного уравнения

$\frac{tg\mu L}{\mu }=\frac{b_1+b_2}{{\mu }^2-b_1{b}_2}$. (6)

Тепловое поле пластины в изображении по Лапласу по временной координате определяется следующим образом:

$\begin{array}{l}Q(x,p)=\underset{0}{\overset{L}{\int }}W(x,\xi ,p)\omega (\xi ,p)d\xi =\\ =a_C\left[b_2{Q}_2\left(p\right)W(x,L,p)-b_1{Q}_1\left(p\right)W\left(x\mathrm{,0,}p\right)\right]+\underset{0}{\overset{L}{\int }}W(x,\xi ,p)Q_0\left(\xi \right)d\xi .\end{array}$ (7)

Выражение (7) позволяет осуществить расчет температуры в произвольной точке x пластины, в зависимости от поведения температур сред на границах. В случае переменных коэффициентов аналитического решения задачи (1) – (3) не существует. Далее предлагается подход к использованию аналитического решения (7) для реализации численного расчёта температуры пластины при изменяющихся во времени коэффициентах теплопередачи.

Решение задачи при переменных коэффициентах

Выражение (5) может быть представлено в виде

$W(x,\xi ,p)=\sum _k\frac{K_k(x,\xi )}{p+a_C{\mu }_k^2}$, $k=\mathrm{1,2,...}$, (8)

где

$K_k(x,\xi )=\frac{{\varphi }^{*}({\mu }_k,x){\varphi }^{*}({\mu }_k,\xi )}{{||{\varphi }^{*}({\mu }_k,x)||}^2}$.

С учётом (8) решение (7), при $Q_0\left(x\right)=0$, может быть представлено в виде бесконечной суммы параллельных апериодических звеньев первого порядка:

$Q(x,p)=\sum _k{Q}_k(x,p)=\sum _k\left[b_2{K}_k(x,L)Q_2\left(p\right)-b_1{K}_k\left(x\mathrm{,0}\right)Q_1\left(p\right)\right]\frac{1}{p+a_C{\mu }_k^2}$. (9)

Каждое слагаемое $Q_k(x,p)$ может быть представлено в виде структурной схемы (рис. 1).

 

Рис. 1. Структурная схема k-го слагаемого решения с постоянными коэффициентами

 

Изменение параметров задачи (1) – (3) в произвольный момент времени t может рассматриваться как начало нового динамического процесса, с новыми параметрами и начальным температурным распределением, существовавшем в момент изменения [19]. Структура, представленная на рис. 1 иллюстрирует, что температура в точке x в произвольный момент времени будет определяться суммой температур на выходах интеграторов. Таким образом, при реализации расчёта в пакетах моделирования динамических систем, текущее состояние будет использоваться как начальное при новых значениях коэффициентов. Необходимо только реализовать вычисление коэффициентов модели в зависимости от изменившихся параметров краевой задачи.

Изменение коэффициентов теплопередачи влияет на значения корней трансцендентного уравнения (6). После ввода обозначений

$f_{1K}(\mu ,x)=\frac{{\varphi }^{*}(\mu ,x){\varphi }^{*}\left(\mu \mathrm{,0}\right)}{{||{\varphi }^{*}\left(\mu \right)||}^2}$,

$f_{2K}(\mu ,x)=\frac{{\varphi }^{*}(\mu ,x){\varphi }^{*}(\mu ,L)}{{||{\varphi }^{*}\left(\mu \right)||}^2}$,

структурная схему для слагаемого $Q_k(x,t)$ преобразована следующим образом (рис. 2). Сигналы b₁(t), b₂(t) определяются так:

$b_1\left(t\right)={\alpha }_1\left(t\right)/\lambda \left(t\right)$,

$b_2\left(t\right)={\alpha }_2\left(t\right)/\lambda \left(t\right)$.

 

Рис. 2. Структурная схема k-го слагаемого решения с переменными коэффициентами

 

Семейство функций ${\mu }_k(b_1,b_2)$, $k=\mathrm{1,2,...}$, реализует вычисление корней трансцендентного уравнения (6).

Для реализации в компьютерных пакетах моделирования динамических систем более удобным может быть представление модели в пространстве состояний вектора $\mathbf{\theta }={\left[{\theta }_i\right]}_{N\times 1}$:

$\left\{\begin{array}{l}\frac{d\mathbf{\theta }\left(t\right)}{dt}=\mathbf{A}\left(b_1\right(t),b_2(t),t)\mathbf{\theta }\left(t\right)+\mathbf{B}\left(b_1\right(t),b_2(t),x)\mathbf{u}\left(t\right);\\ Q(x,t)=\mathbf{C\theta }\left(t\right).\end{array}\right.$ (10)

Здесь

$\mathbf{A}(b_1,b_2,t)=-\mathrm{diag}\left(\left[\begin{array}{cccc}a\left(t\right){\mu }_1^2(b_1,b_2)& a\left(t\right){\mu }_2^2(b_1,b_2)& \mathrm{...}& a\left(t\right){\mu }_N^2(b_1,b_2)\end{array}\right]\right)$,

$\mathbf{B}(b_1,b_2,x)=\left[\begin{array}{cc}\frac{f_{1K}\left({\mu }_1\right(b_1,b_2),x)}{f_T\left({\mu }_1\right(b_1,b_2\left)\right)}& -\frac{f_{2K}\left({\mu }_1\right(b_1,b_2),x)}{f_T\left({\mu }_1\right(b_1,b_2\left)\right)}\\ \frac{f_{1K}\left({\mu }_2\right(b_1,b_2),x)}{f_T\left({\mu }_2\right(b_1,b_2\left)\right)}& -\frac{f_{2K}\left({\mu }_2\right(b_1,b_2),x)}{f_T\left({\mu }_2\right(b_1,b_2\left)\right)}\\ \mathrm{...}& \mathrm{...}\\ \frac{f_{1K}\left({\mu }_N\right(b_1,b_2),x)}{f_T\left({\mu }_N\right(b_1,b_2\left)\right)}& -\frac{f_{2K}\left({\mu }_N\right(b_1,b_2),x)}{f_T\left({\mu }_N\right(b_1,b_2\left)\right)}\end{array}\right]$,

$\mathbf{C}={\left[\begin{array}{cccc}1& 1& \mathrm{...}& 1\end{array}\right]}_{1\times N}$,

$\mathbf{u}\left(t\right)={\left[\begin{array}{cc}{Q}_1\left(t\right)& Q_2\left(t\right)\end{array}\right]}^{\mathrm{T}}$,

N – количество слагаемых бесконечного ряда (9), учитываемых в численной реализации модели.

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

Описанный подход применён при реализации модели температуры стенки барабана парового котла ТГМ-84 (Е420/140ГМ ТКЗ). Ранее [20] была проведена идентификация коэффициентов теплообмена, принятых постоянными в течение всего времени нагрева. При этом получены удовлетворительные результаты поведения модели на этапе растопки, однако, при выходе на рабочий режим наблюдается существенное отклонение модельной температуры от реальной.

Поиск корней трансцендентного уравнения (6) реализован в пакете Matlab^® с помощью функции fminbnd, которая осуществляет поиск минимума функции одной переменной на заданном интервале. Исходя из результатов идентификации [20], предположено, что коэффициент ${\alpha }_1\in \left[0.01;100\right]$ Вт/(м²К), ${\alpha }_2\in \left[10;200\right]$ Вт/(м²К), $\lambda =48$ Вт/(м·К), L=0.09 м. На каждом отрезке взято по 20 опорных точек, для ${\alpha }_1$ – равноотстоящих логарифмически, для ${\alpha }_2$ – линейно.

Уравнение (6) имеет особое решение, в случае, когда ${\mu }^2=b_1{b}_2$, что должно быть учтено в алгоритме формирования интервала отыскания корня. Условия для определения интервалов могут быть сформулированы следующим образом:

1. если $b_1{b}_2\le \pi /\left(2L\right)$, то первый корень ${\mu }_1$ располагается на отрезке $\left[b_1{b}_2;\pi /\left(2L\right)\right]$, последующие корни располагаются на интервалах $\left(\pi \left(2i-1\right)/\left(2L\right);\pi \left(2i+1\right)/\left(2L\right)\right)$, $i=\mathrm{2,3,...}$;
2. если $b_1{b}_2>\pi /\left(2L\right)$, то на интервале $\left(0;\pi /\left(2L\right)\right]$ корня нет, расположение корней подчиняется следующему закону: на интервалах, для которых выполняется условие

$C_1=b_1{b}_2\notin \left(\pi \left(2i-1\right)/\left(2L\right);\pi \left(2i+1\right)/\left(2L\right)\right]$, $i=\mathrm{2,3,...}$,

располагается один корень; на интервале j, для которого условие C₁ не выполняется, располагается два корня, первый из которых лежит в интервале $\left(\pi \left(2j-1\right)/\left(2L\right);b_1{b}_2\right)$, а второй – на отрезке $\left[b_1{b}_2;\pi \left(2j+1\right)/\left(2L\right)\right]$.

На рис. 3 представлены поверхности отражающие зависимость ${\mu }_k(b_1,b_2)$ для первых трёх корней уравнения.

 

Рис. 3. Зависимости первых трёх корней уравнения (6) от b₁, b₂

 

По результатам анализа вида зависимостей ${\mu }_k(b_1,b_2)$, было принято решение при реализации модели для первого корня применить линейную интерполяцию по рассчитанным данным, для остальных – использовать билинейную аппроксимацию.

Модель (10) была реализована в пакете Matlab^® Simulink^® с помощью S‑функции (S-Function) – специального механизма пакета, позволяющего описать произвольный алгоритм расчёта модели объекта, представленного в пространстве состояний. В алгоритме S-функции реализован расчёт зависимостей ${\alpha }_1\left(T_{wall}^M\right)$, ${\alpha }_2\left(T_{wall}^M\right)$, $\lambda \left(T_{wall}^M\right)$, $с\left(T_{wall}^M\right)$, где $T_{wall}^M$ – рассчитанная (модельная) температура стенки в точке контроля $x=0$, $T_{wall}^M=Q\left(\mathrm{0,}t\right)$. Зависимости теплоёмкости и теплопроводности аппроксимировались кусочно-линейной зависимостью по табличным данным для низколегированной стали 22К [21]. Зависимости коэффициентов теплопередачи определены в процессе идентификации модели. Координаты точки $x=0$ соответствуют внешней теплоизолированной границе стенки барабана котла, $x=L$ – внутренней границе стенки, контактирующей с пароводяной смесью.

Идентификация модели

Идентификация модели проводилась по реальным данным, полученным в процессе растопки котла ТГМ-84. В качестве исходных данных использовались: температура воды в барабане котла $T_B\left(t\right)$, температура стенки барабана котла $T_{wall}\left(t\right)$. На первом этапе коэффициенты теплопередачи были определены в виде линейной зависимости от температуры стенки:

${\alpha }_1\left(T_{wall}^M\right)={\alpha }_{11}{T}_{wall}^M+{\alpha }_{10}$,

${\alpha }_2\left(T_{wall}^M\right)={\alpha }_{21}{T}_{wall}^M+{\alpha }_{20}$.

Подбор коэффициентов осуществлялся путём минимизации функционала

$J\left(\Delta \right)=\underset{0}{\overset{t_1}{\int }}\left|T_{wall}^M(\Delta ,t)-T_{wall}\left(t\right)\right|dt\to \underset{\Delta }{\min }$,

где $\Delta =\left[{\alpha }_{10},{\alpha }_{11},{\alpha }_{20},{\alpha }_{21}\right]$ – вектор варьируемых параметров. Область определения для каждого параметра определялась таким образом, чтобы значения ${\alpha }_1\left(T_{wall}^M\right)/\lambda \left(T_{wall}^M\right)$, ${\alpha }_2\left(T_{wall}^M\right)/\lambda \left(T_{wall}^M\right)$ не выходили за пределы областей определения для b₁ и b₂, диапазон изменения температуры стенки взят из экспериментальных данных, $T_{wall}^M\in \left[\mathrm{20,}320\right]$ °С.

В результате решения задачи получены следующие результаты:

1. коэффициент теплопередачи на внешней стороне стенки может быть принят постоянным ${\alpha }_1\left(T_{wall}^M\right)=10.5$ Вт/(м²К);
2. линейная зависимость для коэффициента теплопередачи на внутренней стороне стенки не обеспечивает требуемого качества поведения модели: на первом временном интервале процесса растопки, до достижения температуры металла 150 °C, модельная температура превышает фактическую, на втором интервале – отстаёт;
3. при выходе на установившийся режим модельная температура становится существенно выше фактической.

На втором этапе идентификации было принято решение аппроксимировать ${\alpha }_2\left(T_{wall}^M\right)$ кусочно-линейной зависимостью, заданной опорными точками. Для повышения скорости расчёта, идентификация проводилась в два этапа: на первом и на втором временных интервалах. В результате идентификации были найдены оптимальные по критерию минимума интеграла модуля ошибки значения коэффициентов теплопередачи в узловых точках (рис. 4). Сшивка интервалов производилась по значению вектора состояний в конце первого интервала.

 

Рис. 4. Зависимость коэффициента теплопередачи от температуры стенки барабана

 

Полученная зависимость коэффициента теплопередачи от температуры обеспечивает высокую точность модели, но совершенно нефизична: модель переобучена. Этому есть несколько объяснений: 1) нелинейное поведение датчиков температуры воды в барабане котла, а также стенки барабана; 2) влияние на коэффициент теплопередачи дополнительных факторов. Действительно, в процессе растопки машинист имеет возможность воздействовать на температуру стенки барабана с помощью продувки котла, подачи пара в линию обогрева. Кроме того, очевидно, что коэффициент теплопередачи существенно уменьшается при выходе котла на рабочий режим. Было сделано предположение, что коэффициент ${\alpha }_2$ зависит не только от температуры пароводяной смеси в барабане, но и от $F_{SW}$ – расхода питательной воды в барабан котла. Поэтому на третьем этапе идентификации было принято решение уменьшить число опорных точек (рис. 4), для зависимости коэффициента теплопередачи от температуры ${\alpha }_2\left(T_{wall}^M\right)$ и добавить влияние расхода питательной воды в виде:

${\alpha }_2\left(T_{wall}^M,F_{SW}\right)={\alpha }_2\left(T_{wall}^M\right)-K_F\max \left(\mathrm{0,}{F}_{SW}-F_0\right)$.

Это позволило учесть пропорциональное уменьшение коэффициента теплообмена при увеличении расхода питательной воды выше $F_0=100$ т/час, $K_F=0.58$. Результаты моделирования приведены на рис. 5 – 7. Точность модели по интегральному критерию модуля ошибки возросла в 2.8 раза по сравнению с моделью с постоянным коэффициентом [20].

 

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

 

Рис. 6. Рассогласование фактической и модельной температур

 

Рис. 7. Изменение расхода питательной воды и коэффициента теплопередачи в процессе моделирования

 

Заключение

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

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

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

Об авторах

Иван Александрович Данилушкин

Самарский государственный технический университет

Автор, ответственный за переписку.
Email: idanilushkin@mail.ru
ORCID iD: 0000-0002-2009-2948
SPIN-код: 7289-5030
Scopus Author ID: 57202161857
ResearcherId: D-6290-2014

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

Россия, 443100, г. Самара, ул. Молодогвардейская, 244

Сергей Александрович Колпащиков

Самарский государственный технический университет

Email: skolpaschikov@mail.ru

кандидат технических наук, заведующий кафедрой автоматики и управления в технических системах

Россия, 443100, г. Самара, ул. Молодогвардейская, 244

Илья Сергеевич Левин

Самарский государственный технический университет

Email: Levin.is@samgtu.ru

(Ph.D. (Techn.)), Associate Professor

Россия, 443100, г. Самара, ул. Молодогвардейская, 244

Список литературы

Дополнительные файлы