Идентификация параметров внешнего теплообмена в зоне вторичного охлаждения МНЛЗ Текст научной статьи по специальности «Математика»

УДК 681.5 : 51-74

В.Н. ТКАЧЕНКО, д-р техн. наук,

А.А. ИВАНОВА,

Г.Р. ВАСИЛЯН

ИДЕНТИФИКАЦИЯ ПАРАМЕТРОВ ВНЕШНЕГО ТЕПЛООБМЕНА

В ЗОНЕ ВТОРИЧНОГО ОХЛАЖДЕНИЯ МНЛЗ

В роботі розглянуто задачу ідентифікації коефіцієнту тепловіддачі у вторинній зоні охолодження машини безперервного лиття заготовок. Нестаціонарні процеси теплообміну описуються нелінійними параболічними рівняннями у частинних похідних. Положення невідомої межі розподілу фаз визначається умовами Стефана. Запропоновано метод розв’язання поставленої задачі за допомогою метода найменших квадратів. Також розроблено алгоритм розв’язання цієї задачі на основі методу стохастичної апроксимації. Наведені результаті розрахунків за обома цими методами і виконано їх аналіз.

The identification problem of heat-transfer coefficient in the secondary cooling zone of continuous casting machine is considered. Non-stationary processes of internal heat exchange are described by the nonlinear parabolic partial differential equations. The solution of the problem using the least squares method is proposed. The solution algorithm based on stochastic approximation method is developed. The analysis of numerical results is performed.

Введение. При построении математических моделей технологических процессов возникают трудности в подборе некоторых параметров процесса. В частности при моделировании процесса тепломассопереноса в слитке машины непрерывного литья заготовок (МНЛЗ) каждый раз остаётся открытым вопрос выбора коэффициента теплоотдачи (КТ) на поверхности слитка в зоне вторичного охлаждения (ЗВО). Связано это с тем, что на величину КТ влияет множество различных факторов. Кроме того, значение КТ может достаточно сильно изменяться во времени и по пространственным координатам. Таким образом, возникает задача идентификации распределённого параметра.

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

Постановка задачи. В системе координат, привязанной к неподвижной конструкции криволинейной МНЛЗ, рассматривается тепловое поле движущегося стального слитка и стенок кристаллизатора. Общая схема объекта представлена на рис. 1.

На участке слитка внутри кристаллизатора задано уравнение тепломассопереноса и граничные условия [1]. На границе раздела фаз заданы условия равенства температур и условия Стефана, а также граничное и начальное условия для границы фазового перехода.

Рис. 1. Схема машины непрерывного литья заготовок (МНЛЗ). 1 - жидкая фаза; 2 -

твёрдая фаза; 3 - точка окончательного затвердевания слитка; 4 - форсунки, распыляющие водо-воздушную смесь; 5 - обжимные и приводные ролики

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

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

Уравнение тепломассопереноса для слитка на криволинейных участках МНЛЗ выглядит следующим образом:

1

+в „ о =----------------------

Эх Эф с(Т, г, ф)р(Т, г, ф)

д

Эг

ЭТ Л 1 Э

— I ЦТ,г,ф)— І + —• —I ЦТ,г,ф)— 1 +

Эг )

Эф

ЭТ Л Х(Т, г, ф) ЭТ

Эф І)

Эг

(1)

где 9т - угловая скорость движения слитка на т-м криволинейном участке. Условия для неизвестной границы на криволинейных участках:

Т (Т’Г’ Ф)| г=?12_ (т,„) = Т (Т’Г’ ^1 г

'■г=6 2+ (г,ф) Ткг ;

\дТ КТ, г, ф) —

ЭН

\дТ

-Х(Т, г, ф)—

ЭН

= ЦРкг I 0т (г)+ ^|;

г

г

+

Si(0, ф) = ^і0(ф);

\(Т лдТ

Х(Т, r, ф)— = —

(2)

£ 2 (0, Ф) = £ 20 (Ф),

где |1(ф) и 12(ф) - границы раздела фаз.

Граничные условия на криволинейных участках:

где ai(Gm(T), ф), Ote(Gm(r), ф), Cm, Свш, Ты, Твш - коэффициенты теплоотдачи на поверхности слитка, приведённые коэффициенты излучения от поверхности слитка, температура окружающей среды в ш-й секции ЗВО по внутреннему и по внешнему радиусам соответственно, Gm(z) - расход воды в m-й секции.

Аналогично задаётся уравнение тепломассопереноса на прямолинейном участке и граничные условия на прямолинейном участке. Если жидкая фаза продолжается дальше точки выпрямления, то на прямолинейном участке ЗВО также задаются условия для неизвестной границы раздела фаз.

Считаем, что в конце прямолинейного участка тепловой поток равен нулю:

Причём начальная температура есть непрерывная функция на всей области слитка и на стенке кристаллизатора.

Требуется по имеющейся информации о температуре на поверхности слитка определить коэффициент теплоотдачи аі(Ош(т), ф), ал(От(т), ф) в ЗВО

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

(5)

(6)

МНЛЗ.

условие, что легко можно проверить, используя для решения поставленной задачи метод прямого обращения [2].

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

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

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

Рассмотрим участок слитка внутри первой секции ЗВО. Коэффициент а.1(От(т), ф) имеет специальное распределение вдоль поверхности слитка. Известно, что на участке, накрываемом факелом форсунки, его можно приблизить параболической функцией, которая приобретает максимальное значение в точке, соответствующей координате сопла форсунки, а на остальных участках - константой. В одной секции установлены однотипные форсунки, которые дают одинаковый водо-воздушный факел, следовательно, и коэффициент теплоотдачи - одна и та же парабола, сдвинутая вдоль оси абсцисс (рис. 2).

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

было сказано, а = const, а на участках, подвергающихся принудительному охлаждению, будем искать а(ф) в виде

Р 2

а(ф) = ас Р Ф +а р ■

W

(7)

Рис. 2. Распределение коэффициента теплоотдачи вдоль поверхности слитка внутри

одной секции

Рассмотрим сначала участки, на которых а(ф) = ас = const. Обозначим множество узлов ф,-, в которых КТ считается постоянным, К. Множество остальных узлов, где КТ распределяется согласно параболическому закону, обозначим В.

Выберем следующую конечно-разностную аппроксимацию граничного условия (3):

X

T,2 - 4Тд +

7,0

2q

= аc(Th -Tt,0) + Ch(T4 -Тф),

(8)

где q - шаг конечно-разностной сетки по радиусу п.

Отсюда получаем формулу для невязки тепловых потоков на границе:

А = X

T-2 - 4T,1 + 3T,0

t,0

2q

-ChTx -T40)-ac(Th -T,0).

Обозначим

T. 2 — 4Tj 1 + 3Tj 0 4 4

Р =Ь i,0 —---------fg------ — Ch — Tj4 ), P2i = Th — T,0.

Требуется найти величину a^, чтобы сумма квадратов невязок была минимальной, т.е. чтобы выполнялось условие

S = £(Ри — acP2i)2 ^ min, Vi : 9i е К .

Необходимое условие существования экстремума S(o(c):

_с$_

Са„

= -2£ Р2і Р-а сР2і) = 0.

Отсюда находим ас

Каждому узлу ф,- из множества В поставим в соответствие точку у,■ на отрезке [- V, w] таким образом, чтобы |у,| равнялось расстоянию от соответствующего ф,- до координаты ближайшей форсунки р,. Из (7) и (8) получаем невязку

А = Х

Ті,2 - Т + 37!,о

і,0 '

2д

-Сл(Т74 -Ті4о)-

а,

- (а с -^ У? +а р )(Т7і - Т,о). м>

Аналогичным способом находим ар

‘с XР2

-X Р1 і Р2

2

*Г -1

V У

X Р2

Следует также отметить, что в определённых при помощи МНК значениях а и ар выполняется достаточное условие существования минимума функции & Легко проверить, что частные производные второго порядка £ по каждому их этих параметров строго больше нуля.

На рис. 3 представлены сравнительные результаты вычислений. В качестве теплофизических параметров для модели были выбраны данные процесса непрерывной разливки стали (для марки ст 40), ширины сляба 1 м, полутолщины сляба равны 0,1 м и скорости движения слитка V =1 м/мин. Здесь хорошо видно, что решение, полученное методом прямого обращения, является неустойчивым и непригодным для практического использования. Вторая кривая представляет сплайн аппроксимацию, которая является результатом решения той же задачи методом наименьших квадратов.

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

а. =

с

2

а„ =

р

2

2

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

Аналогично определяется коэффициент теплоотдачи для остальных секций зоны вторичного охлаждения.

Оперативная настройка коэффициента теплоотдачи.

После начальной настройки необходимо обеспечить адаптацию параметров модели во время обычной работы МНЛЗ, когда информация о тепловом состоянии ограничивается показаниями температуры иногда всего только в одной точке на поверхности слитка. Кроме всего, подстройка параметра должна выполняться в реальном масштабе времени. Такие алгоритмы могут быть основаны на методе стохастической аппроксимации [12].

Длина слитка, м

Рис. 3. Результаты определения коэффициента теплоотдачи: 1 - методом прямого обращения; 2 - методом наименьших квадратов

Температура на поверхности слитка измеряется через равные небольшие промежутки времени. Обозначим измеряемую температуру Т/. Во время работы МНЛЗ производственный процесс разливки параллельно моделируется с помощью описанной выше математической модели в управляющей вычислительной машине. Температуру в соответствующей точке, рассчитанную по модели, обозначим 7}. Необходимо по отклонениям расчетной температуры от измеренной корректировать параметры модели так, чтобы свести к минимуму эти отклонения. Трудность решения данной задачи состоит в том, что измерения температуры искажены помехой.

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

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

Алгоритм настройки параметра имеет вид:

______

а }+1 =а} - к} 7 - Т} ),

где а} - }-е приближение коэффициента теплоотдачи а, к} - специальная последовательность чисел, которая удовлетворяет следующим условиям:

¡1Ш к} = 0 , V к} =ж, V (ki) } } ¿.л

}=1 }=1

< .

Таковой является, например простейшая последовательность

а

к} =- ,

1 Ь + }

где а, Ь е Я, а > 0. Подбирая числа а и Ь, а также другие последовательности, удовлетворяющие вышеуказанным условиям, можно изменять скорость сходимости алгоритма. В [12], например, рекомендуется поддерживать к} постоянным до изменения знака невязки (7}* - Т), изменяя затем к} так, чтобы удовлетворить вышеупомянутым ограничениям.

Критерием точности подстройки параметра служит вхождение последних полученных N приближений в окрестность истинного значения КТ. В качестве истинных значений КТ служат величины, которые были подобраны экспериментально при решении прямой задачи - моделирования теплового поля МНЛЗ [1].

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

Примеры. Рассмотрим простую последовательность чисел

к} = -, } = 1,2,3,...

для различных значений коэффициента а. На рис. 4 представлены результаты работы этого алгоритма. При значениях коэффициента а <1 наблюдается слишком медленная сходимость. Если выбрать а = 1, то приблизительно после 150-й итерации значение подстраиваемого параметра будет находиться в достаточно малой окрестности истинного значения. При а = 2 траектория подстройки параметра отражает колебания с затухающими амплитудой и частотой, и уже после 80-ти итераций параметр настроен. При увеличении а > 2 колебания становятся слишком большими и алгоритм, хотя и достаточно быстро находит истинное значение параметра, но затем "проскакивает" его. В этом случае также наблюдаются колебания с затухающими амплитудой и частотой, но для подстройки требуется значительно большее число итераций.

Число тактов

Рис. 4. Траектории подстройки параметров, характеризующие отклонение подстраиваемого параметра от истинного значения.

1 - а = 0,5; 2 - а = 1; 3 - а = 2; 4 - а = 3.

Таким образом, можно сделать вывод, что для выбранной

последовательности наилучшими являются значения коэффициента а из отрезка [1, 2].

Исследуем теперь влияние величины Ь на скорость сходимости. На рис. 5 показаны траектории подстройки параметра при различных значениях Ь. Видно, что увеличение Ь не имеет смысла в данном случае, так как только замедляет сходимость алгоритма.

Число тактов

Рис. 5. Траектории подстройки параметров, характеризующие отклонение подстраиваемого параметра от истинного значения:

1 - b = 0; 2 - b = 5; 3 - b = 10.

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

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

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

Список литературы: 1. Ткаченко В.Н., Иванова А.А. Анализ температурных полей

криволинейной МНЛЗ на основе математического моделирования / Матеріали 3-ї міжнародної науково-практичної конференції "Прогресивні технології у металургії сталі: ХХІ сторіччя". -Донецьк: ДонНТУ. - 2007. - С. 242-249. 2. Ткаченко В.Н. Моделирование тепловых процессов в автоматизированных системах обработки информации // Вісник Донецького національного університету. Серія А. Природничі науки. - 2002. - № 2. - С. 379-383. 3. Алифанов О.М., Артюхин В.Д., Румянцев С.В. Экстремальные методы решения некорректных задач. - М.: Наука, 1988. - 286 с. 4. КоздобаЛ.А., Круковский П.Г. Методы решения обратных задач теплопереноса. -К.: Наукова думка, 1982. - 360 с. 5.Мацевитый ЮМ. Обратные задачи теплопроводности. В 2-х т. - НАН Украины, Институт проблем машиностроения. - К.: Наукова думка, 2003.

6.МикитенкоН.И. Сопряженные и обратные задачи тепло-массопереноса. - К.: Наукова думка, 1988. - 240 с. 7. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. - М.: Наука, 1974. - 224 с.

Поступила в редакцию 20.09.2007