Автоматизированное проектирование оптимальной теплообменной системы на основе решения задачи процесса теплопереноса Текст научной статьи по специальности «Математика»

4. Книга, Е.В. Принципы построения комбинированной топологии сети для перспективных бортовых вычислительных систем / Е.В. Книга, И.О. Жаринов // Научно-технический вестник информационных технологий, механики и оптики. - 2013. - № 6. - С. 92 - 98.

5. Копорский, Н.С. Система бортовой картографической информации пилотируемых летательных аппаратов. Основные принципы построения / Н.С. Копорский, Б.В. Видин, И.О. Жаринов // Сборник трудов Х Международной конференции «Теория и технология программирования и защиты информации». - СПб., 2006. - С. 18 - 23.

6. Костишин, М.О. Оценка точности визуализации местоположения объекта в геоинформационных системах и системах индикации навигационных комплексов пилотируемых летательных аппаратов / [М.О. Костишин и др.] // Научно-технический вестник информационных технологий, механики и оптики. - 2014. - № 1. - С. 130 - 137.

7. Парамонов, П.П. Интегрированные бортовые вычислительные системы: обзор современного состояния и анализ перспектив развития в авиационном приборостроении / П.П. Парамонов, И.О. Жаринов // Научно-технический вестник информационных технологий, механики и оптики. - 2013. - № 2 (84). - С. 1 - 17.

8. Парамонов, П.П. Принцип формирования и отображения массива геоинформационных данных на экран средств бортовой индикации / [П.П. Парамонов и др.] // Научно-технический вестник информационных технологий, механики и оптики. - 2013. - № 6. - С. 136 - 142.

9. Парамонов, П.П. Реализация структуры данных, используемых при формировании индикационного кадра в бортовых системах картографической информации / [П.П. Парамонов и др.] // Научно-технический вестник информационных технологий, механики и оптики. - 2013. - № 2 (84). - С. 165 - 167.

10. Парамонов, П.П. Структурный анализ и синтез графических изображений на экранах современных средств бортовой индикации на плоских жидкокристаллических панелях / [П.П. Парамонов и др.] // Авиакосмическое приборостроение. - 2004. - № 5. - С. 50 - 57.

11. Парамонов, П.П. Теория и практика статистического анализа картографических изображений в системах навигации пилотируемых летательных аппаратов / П. П. Парамонов, Ю.А. Ильченко, И.О. Жаринов // Датчики и системы. - 2001. - № 8. - С. 15 - 19.

УДК 681.5.015:66.021.4:678.074

С.Ю. Осипов, Ю.Р. Осипов, О.А. Панфилова, И.А. Новожилов

АВТОМАТИЗИРОВАННОЕ ПРОЕКТИРОВАНИЕ ОПТИМАЛЬНОЙ ТЕПЛООБМЕННОЙ СИСТЕМЫ НА ОСНОВЕ РЕШЕНИЯ ЗАДАЧИ ПРОЦЕССА ТЕПЛОПЕРЕНОСА

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

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

The article presents the solution of boundary value problem of heat transfer in heat in the object of arbitrary form by the finite element method. The possibility of optimization of thermal modes vulcanizing equipment by creating a control system of temperature control is investigated. The mathematical model of the object with distributed parameters for predicting thermal processes at vulcanization of rubber products is developed.

Finite element method, temperature, thermal regime, mathematical model, heat treatment, functional, object with distributed parameters, loss of production, automated control system.

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

_дГХдГ\ + _д. fx— 1+ - 0 (1)

дх У дх J ду У ду J v

с нулевыми начальными условиями: T (x, y,0)- 0,

idT + T

X--+ al

dn

(2)

при граничных условиях третьего рода:

- qs, (3)

где X - теплопроводность; Т - температура; а - коэффициент теплоотдачи на границе; , д5 - объ-

L

емная и поверхностная плотности мощности источников теплоты.

Величины X, а, qv , qS могут быть заданы в виде произвольных функций координат х, у, в том числе и кусочно-непрерывных функций.

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

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

функционал I [Т (х, у)] вида

где F¡(п), Fjn, Fk[n| - линейные функции координат x, у, равные единице в узлах i, ] или k соответственно и равные нулю в двух других узлах. Индекс (п) указывает, что функции формы относятся к п-му элементу.

Таким образом, для функции формы Fг■n\x,у) должны выполняться равенства

F(п)( х,, у, ) = 1, F(п)( xJ, у ) = F(п)(хк, Ук ) = 0. (5)

Аналогичные соотношения имеют место и для двух других функций Fjn (х, у), Fkn (х, у). Используя условия (5), можно выразить коэффициенты а^т , , Ст (т = ,, J, к) через координаты узлов ,, J, к и определить функции формы. Из условий (5) следует:

Лп)

' [Т (x, У )] = Ц

<]2 + >{!)2 -2,Т

йхйу +

а, =(х1Уь - хьУ1) /2 Ъ = (( -Ук)/2S ;

С =(хк - XJ)/ 2

I aJ =(( - хУ)/ 2 Щ = (Ук - У, )/2S ; к = (х, - хк )/2

■ДаТ2 - 2qsT).

(4)

В методе конечных элементов приближение для искомой функции Т(х, у) отыскивается в виде

Т (^ У )~£ ат/т (X, У)

где ат - неизвестные постоянные коэффициенты, а /т (х, у) - известные функции пространственных координат.

Разобьем гуммированный объект на N треугольных элементов и введем М узлов во всех вершинах треугольников. Причем вершины одних треугольников не должны находиться на сторонах других треугольников, т. е. каждый элемент должен содержать три узла. Присвоим сквозную нумерацию всем элементам (п = 1, ..., N и всем узлам (т = 1, ..., М). Номер узла не связан с номером элемента.

Искомыми величинами в методе конечных элементов являются приближенные значения температуры tт в узлах т = 1, ..., М. Согласно основной идее метода, распределение температуры в каждом элементе представляется в виде суммы, в которую входят три функции формы элемента, умноженные на приближенные значения температуры в его трех узловых точках. Поэтому распределение температуры в п-м треугольном элементе t(п)(х, у) имеет вид:

= (хгУл - х1Уг) ^

К = (у - у ) ,

- X, )/2S

=

где S - площадь треугольника, значение которой может быть вычислено через координаты узлов:

s = (х1Уь - хьУ1 + Ух - Укх, + хкУ<- х1Уг)/2.

Градиент температуры в каждом элементе имеет постоянное значение, и производные по координатам определяются соотношениями:

дt

Ы

— = + bJtJ + , -т- = + cJtJ + СА.

дх ду

Система уравнений для определения температур в узлах } ММ=1 составляется на основании условий

минимума функционала (4). Этот функционал можно представить в виде суммы интегралов по всем элементам:

I = 11

(п)

(6)

^ = |

+х(| I - ^

йхёу + | (аt2 - 2qst) <Н .

¿п)

(7)

¿п) (х,у) = tlFl(п) (х,у) + tJ (х, у) + п) (х, у)

I = * р(п)

(п).

г (п)/

Второй интеграл в (7) вычисляется лишь для граничных элементов вдоль тех сторон, которые прилегают к внешней границе Ь области Б. Условия ми-

s

т=1

п=1

нимума функционала можно с учетом (6) переписать в виде

д

N д1 (">

-II1 = т = 1, М (8)

И=1 д(

Вычислим сначала интеграл по площади элемен-

та:

д_

I )' -

йхйу =

= | 2+ Ь/, + Ь^) + + х ( + с/, + Скгк)с - ч^ (х у)] ёхёу = = 25 ^(Ь^, + ЬЬ^ + ЬЪЛ )-

- х +с с+ ссЛ)- ч/3].

Если стороны элемента принадлежат границе области, то в функционале следует учесть интеграл по этим сторонам. Пусть граничной является сторона Ь, между узлами ' и у. Тогда

д

(^2 -) =

Ь, '

= | 2 [а(Ч + )- д,^ ] Л = 2Ь

3 6 2

Л

Окончательно получаем для производной от функционала I(п) по температуре ч, следующее выражение:

дГ

(п)

д^

= Х5(п) (ч + ЬЬ/, +ЬЬк1к +С2г 1 +ссу1у + С'С^к) +

+ аЬу (ч , /3 + ^ /6) + аЬл (ч ., /3 + 1к /6) - ^ 5(п) /3 -- ЪЬЦ /2- ЪЦк /2.

Причем слагаемые с Ь, и Ь,к записывают лишь в

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

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

ОТ = Ф,

где О - глобальная матрица теплопроводности размером МхМ; Т - вектор-столбец искомых значений температур в М узлах; Ф - глобальный вектор-столбец тепловых потоков.

О = О,.,.

Т =

' Ф ■

, Ф =

)м _ ФМ _

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

1 5

Т (у, т)= — | Т (х, у, т)х, 0 < у < 1, 0 <т<Т ,

5 0

при у = I, было наименьшим в смысле (1). При этом гарантируется достижение заданных свойств в гуммированном изделии. Тем не менее, необходимо создание автоматизированной системы управления технологическим процессом (АСУ ТП) тепловой вулканизации на базе адаптивных математических моделей, работающих в ускоренном и в реальном масштабах времени [3]. Математическое моделирование основано на тождественности в математическом смысле различных по физической природе процессов, протекающих в ОРП и модели, поэтому, изучая тот или иной вариант разрабатываемого режима работы на модели, можно заранее учесть аналогичный режим работы ОРП. Таким же путем, используя моделирование, можно применять математические модели ОРП для прогнозирования тепловых процессов в АСУ ТП вулканизации гуммированных изделий в ускоренном масштабе времени.

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

5

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

Градиентные методы решения задач оптимального управления ОРП, описываемые параболическими уравнениями, рассмотрены в [4], [5]. В соответствии с этими методами требуется, регулируя температуру внешней среды изделия управлением и () , сделать

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

Температура внешней среды для однородного нагреваемого стержня 0 < х < I

«И+1 {t) = un {t) + аn [«. {t) - un {)]

(9)

где а < и () < Ь - числа, выражающие крайние допустимые значения температуры; ип +1 (/) - требуемое управление в п + 1 приближении; ип () - вспомогательное приближение, определяемое как

т т

|а2уф((,t,ип)ип () Ж = тт|а2уф((,t,ип)ип () Ж ,

0 0

где Т, а, и, I - некоторые положительные постоянные коэффициенты, а ф^) - является решением краевой задачи:

d2 У

dy__ __

dt dx2

{x, t)e Q ;

d У dx

= 0, 0 <t <x ;

d У dx

= -уф{/, t), 0 < t < x ; y|t_T = 2[T{x,t,u)-T3 {x)] , 0 < x <l.

Из (3) для множества u = Ju {t): u {t) e L2 [0, x]}.

a < u{t) <b , 0 < t < x следует:

.{t ) =

a, если a 2vT(l, t, un )> 0,

b, если a 2vf{{, t, un )< 0, (10) un {t), если a2vf(l, t, un ) = 0.

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

>{t ) =

a, если a 2vT{{, t, un )>A,

b, если a 2vf{{, t, un )<-A,

un {t), если a2vT{l,t,un)<|a|,

где Д - заданная постоянная величина.

Коэффициент ап > 0, входящий в формулу (9), находится из условия

min qn {a) = qn {an ), qn {a) = J[_un + a{{ _un)] .

После преобразований в соответствии с [3] можно записать

т

Чп (а) = J (ип ) + а| а2у У (I, t, ип ) |\ () - ип (t)] Ж +

0

I 2

+ а21|Т(х,т,ип)-Т(х,т,ип)| ёх . (11)

0

Откуда видно, что чп (а) достигает экстремума в точке

* 0 a„ = —

ja2 v У {{,t,un)[un {t)-un {t)]dt

->0,

2j |T {x, x, un )- T {x, x, un )| dx

это значит, что квадратный трехчлен (11) достигает своего минимума на отрезке 0 < a < 1 в точке

* * ^ л л * л

an = an при an < 1 и в точке an = 1 при an > 1, то есть можно записать an = min JJ ,1j .

Полученное an, подставляется в (9), откуда следует (n + 1)-е приближение для решения задачи. Если an = 0 или T {x, x, йп) = T {x, x, un), квадратный трехчлен (11) вырождается, итерационный процесс прекращается, и управление un {t) = u* будет неполным оптимальным решением задачи.

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

При разработке программного обеспечения для этих систем необходимо формализовать алгоритмы условного прогнозирования режимов работы ОРП применительно к цифровой реализации системы управления. Указанная задача в [3], [5] сформулирована следующим образом. Суммарные потери F {t0, tk) за время работы объекта от начального момента t0 до конечного tk определяются формулой

x=0

tk

F (t0> tk ) = J f (x)dx, (12)

to

где f (т) - текущие потери производства в момент т. Если f (т) можно связать с оперативной информацией об объекте J (т) = {со (т), u (т), v (т)} , Vr e (-«, t), состоящей из полученных к моменту

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

иоп (т) = иоп [J (т)] = arg min F (to, tk).

u [ J (т)]еиоп

С учетом цифровой реализации критерия оптимизации tn = nT (n = 0, 1, 2, ...), где T - фиксированный период. При этом класс допустимых сигналов управления ограничен семейством ступенчатых функций

u (t) = un, Vt e [nT, (n + 1)t] ,

где un - варьируемые переменные.

Оптимальный сигнал un вычисляется как

unü (т) = arg mm {F (и, tn+l)J , unl) + S„ (unl)},

un 1 eu

где Sn (un1) - прогноз будущих потерь в промежутке tn+j, tk. При данной J(tn) и данном un, при условии, что все будущие значения сигнала (7) будут строиться в моменты времени tn+1 =(n + i)T ,

( = 1,2,...,n-k) аналогично (8), т.е. на основании полученной в соответствующие моменты времени информации J (tn + ') будут реализованы оптимальные значения un + t = u°"+i.

Если подынтегральная функция в (12) задана в виде

/ (Г) = / [и ()] + / )] ,

где / и /2 - выпуклые, неотрицательные определенные функции от и и и, то оценить будущие суммарные потери непосредственно в виде функции от условного прогноза можно в виде формулы

_^ (Чп+' , Чп+'+1 )_ =

■ ( , un,2,..., ип'+1 )

гп+,+1

= | [/1 (ип,'+1 )+ /2 (( / Чп, а, ип ))]т.

чп+,

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

Литература

1. Бутковский, А.Г. Методы управления системами с распределенными параметрами / А.Г. Бутковский. - М., 1975.

2. Канторович, Л.В. Приближенные методы высшего анализа / Л.В. Канторович, В.И. Крылов. - М., 1962.

3. Осипов, Ю.Р. Автоматизация технологических процессов гуммировочных производств / Ю.Р. Осипов, С.Ю. Загребин. - М., 2004.

4. Осипов, Ю.Р. Оптимизация распределенной системы управления непрерывным технологическим процессом / Ю.Р. Осипов, А.А. Моисеев, В.В. Павлов // Вестник Оренбургского государственного университета. - 2002. - № 3. -С. 135 - 138.

5. Осипов, Ю.Р. Термообработка и работоспособность покрытий гуммированных объектов / Ю.Р. Осипов. - М., 1992.