сотрудник
Россия
сотрудник
ВАК 05.17.00 Химическая технология
ВАК 05.23.00 Строительство и архитектура
ГРНТИ 20.53 Технические средства обеспечения информационных процессов
Описывается метод нахождения оптимальной плотности источников тепла в движущейся среде внутри области, находящейся в состоянии стационарного теплового баланса с окружающей средой. Получаемые источники имеют минимальную суммарную интенсивность и обеспечивают приемлемое распределение температуры внутри области. Строятся конечномерные аппроксимации задачи, обладающие особым свойством регулярности по функционалу. Это свойство показывает принципиальную возможность численного нахождения квазирешения исходной задачи. Обсуждаются алгоритмы его приближённого нахождения, а также результаты вычислительных экспериментов, проведённых с помощью специально созданного программного комплекса HeatCore.
плотность источников тепла, обратная задача теплопроводности, задача оптимального управления для эллиптических краевых задач, конечномерная аппроксимация, конвекция, тепловой баланс, симплекс-метод
1. Введение. Одним из видов объектов, широко распространённых в различных областях человеческой деятельности, являются системы источников тепла в области пространства, находящейся в состоянии стационарного теплового баланса с окружающей средой. При математическом моделировании таких систем часто возникает связанная с ресурсосберегающими технологиями инженерная задача оптимального распределения источников тепловых полей. В типичной постановке эта задача состоит в таком выборе распределения источников, при котором создаваемое ими температурное поле наименьшим образом отличается от заданного. В качестве критерия оптимизации здесь обычно выступает квадратичный функционал, обладающий свойством коэрцитивности. С математической точки зрения эта задача относится к задачам оптимального управления для эллиптических краевых задач. Существование решений и общие свойства подобных задач для квадратичных целевых функционалов, а также приближённые методы их решения изучались рядом авторов (см. [1,2] и приведённую там библиографию). Эту задачу можно также отнести к обратным задачам теплопроводности (ОЗТ), методы приближённого решения которых рассмотрены в [3]. Однако и здесь речь идет, в основном, о коэрцетивном квадратичном целевом функционале.
На практике при оптимальной организации обогрева жилых и производственных помещений, теплиц и т.д. естественно стремиться получить не заданное температурное поле, а обеспечить некоторый температурный коридор при минимальных энергетических затратах. В настоящей работе и рассматривается задача отыскания такого распределения плотности источников тепла, которое обеспечивает заданный температурный режим при минимальной суммарной мощности источников. При численном решении этой задачи возникает ряд трудностей, и до настоящего времени в полном объёме она практически не рассматривалась. Целевой функционал здесь является линейным и в силу наличия поточечных неравенств, определяющих температурный режим, возникают значительные трудности в установлении существования точного решения. Вообще говоря, точного решения этой задачи может не существовать. В работе [4] введено так называемое квазирешение, которое всегда существует и с прикладной точки зрения вполне приемлемо. В работах [4-6] предлагается способ конечномерной аппроксимации задачи, на основании которого разрабатываются основные алгоритмы, и создаётся методика приближённого нахождения квазирешения. Эти аппроксимации образуют последовательность задач линейного программирования и, как показано в [4,6], эта последовательность обладает особым свойством регулярности по функционалу. Это свойство обеспечивает принципиальную возможность приближённого нахождения квазирешения. Однако построение конечномерных аппроксимаций основано на очень трудоёмкой процедуре определения так называемой обменной матрицы. В настоящей работе изложен один из наименее трудоёмких способов нахождения этой матрицы, и приведены результаты численных экспериментов, показывающих эффективность разработанного метода.
2. Формулировка оптимизационной задачи
Введём функцию , где
- установившаяся температура в точке
ограниченной связной области D, а
- температура окружающей среды. Эта функция должна удовлетворять уравнению
(1)
где c - коэффициент температуропроводности среды, - поле скоростей среды, которое предполагается подчинённым условию
и потенциальным, а
- плотность источников тепла в области D. Для формулировки краевых условий разобьём границу области D на три части:
где
- часть границы, которая является входом среды в область D,
- часть границы, являющейся выходом (стоком) среды, а
- часть непроницаемой границы. Справедливы следующие соотношения:
где - единичный вектор внешней нормали к границе
. Последнее из этих трёх условий означает, что в область D поступает вещество внешней среды с температурой T0. В дальнейшем мы предполагаем, что в нашем процессе присутствует поток тепла через непроницаемую для среды границу, равный
. При этом коэффициент теплопередачи
является функцией, заданной на
. Будем считать выполненными следующие краевые условия:
(2)
Физически эти условия выражают тепловой баланс области D с окружающей средой.
Будем считать, что поле скоростей области D является потенциальным и его потенциал известен. Нас интересует оптимизационная задача нахождения функции
принадлежащей пространству L2(D) и доставляющей минимум линейному функционалу
(3)
при условии, что выполнены равенства (1) и (2), а также неравенства
(4)
Здесь - минимальный и максимальный профили температур, которые считаются непрерывными функциями, задаваемыми в некоторой подобласти
, которая в дальнейшем называется областью контроля температуры. Отметим, что случай
не исключается. В этом случае естественно
Не исключается и случай, когда
. Существование точного решения задачи (3) при условиях (1), (2), (4) в пространстве L2(D) гарантировать нельзя.
Обозначим через нижнюю границу значений функционала
(3), когда
пробегает множество неотрицательных функций из L2(D), удовлетворяющих условиям (1), (2), (4).
Определение 1. Квазирешением оптимизационной задачи (3) при условиях (1), (2), (4) с данным допуском назовём такую функцию
из L2(D), удовлетворяющую ограничениям задачи, для которой выполняется неравенство
Для приближённого нахождения квазирешения дадим более удобную эквивалентную формулировку оптимизационной задачи. Введём новую неизвестную функцию с помощью равенства
Подставляя это выражение в (1), получим следующее уравнение:
(5)
Из условий (2) следует, что функция w должна удовлетворять краевому условию
(6)
где функция задана на границе
области следующим образом:
Оператор , действующий в пространстве L2(D) на достаточно гладкие функции
, подчинённые краевым условиям (6), является самосопряжённым и положительно определённым, а поэтому имеет ограниченный обратный оператор, определенный на L2(D) (см., например, [7, с. 129]). Отсюда следует, что краевая задача (5), (6) имеет единственное решение, которое можно выразить через функцию Грина оператора L. Последнее обстоятельство играет существенную роль в дальнейшем.
Если учесть описанное выше преобразование краевой задачи, то получим следующую формулировку оптимизационной задачи
(7)
где
а G - оператор, обратный по отношению к оператору L. Для задачи (7) тоже можно ставить задачу о нахождении квазирешения , при этом квазирешение исходной задачи
.
3. Конечномерная аппроксимация оптимизационной задачи и её регулярность по функционалу.
Для приближённого нахождения квазирешения построим конечномерную аппроксимацию задачи (7) в виде задачи линейного программирования. Разобьём область D на n частей . Определим подпространство Sn(D)ÌS(D) кусочно-постоянных функций вида
,
. Введём в Sn(D) базис, состоящий из функций
, и
. Тогда
. Разбиение области D мы считаем и разбиением области
, т.е. при некотором натуральном p справедливо равенство
. Рассмотрим на S(D) оператор Gg и обозначим
aij=(mesDi)-1(Gej,ei), ai=(mesDi)-1 (q1,ei), bi=(mesDi)-1 (q2,ei),
где (×,×) - скалярное произведение в L2(D). Заменяя класс функций S(D) подпространством Sn(D), умножая скалярно ограничения на базисные функции , получаем конечномерную аппроксимацию задачи (7)
(8)
Здесь . Самым трудным при построении задачи (8) является нахождение матрицы aij=(mesDi)-1 (Gej,ei), поскольку оператор G явно не задан. Эту матрицу мы и называем обменной матрицей. Её определение равносильно нахождению функций wj = Gej, которые являются решениями уравнений -cDw+(|Ñj|2/(4c))×w=ej при краевых условиях (6). Вопросу о нахождении обменной матрицы в основном и посвящена настоящая работа.
Задачу (8) обозначим через Z0(n). При неограниченном дроблении области D получаем последовательность задач линейного программирования, которую мы и называем конечномерной аппроксимацией задачи (7). Обозначим через (Jn)min минимальное значение целевой функции задачи Z0(n).
Определение 2. Конечномерную аппроксимацию назовем регулярной по функционалу, если справедливо неравенство
,
где g определено перед определением 1.
При наличии регулярности по функционалу решение задачи Z0(n) при достаточно малом можно считать приближённым квазирешением исходной задачи. Действительно, после дискретизации система ограничений означает, что неравенства в исходной задаче (7) практически выполнены, т. к. они удовлетворяются в среднем по Di. При этом значения (Jn)min при достаточно больших n не превосходят g + e.
В работе [6] доказано, что при m ≤ 3 и условиях
1)
2) (9)
конечномерная аппроксимация является регулярной по функционалу.
Назовем разбиение области D, состоящее из подобластей Di, дроблением разбиения, состоящего из подобластей , если это разбиение одновременно является и разбиением каждой области
.
Замечание. Знание обменной матрицы для некоторого разбиения области D позволяет определить обменную матрицу для разбиения, по отношению к которому данное разбиение является дроблением. При этом новая обменная матрица может быть найдена по формуле
(10)
В дальнейшем мы описываем способ определения обменной матрицы (aij) для «тонкого» разбиения области D, а затем с помощью формулы (10) переходим к матрице ( ), отвечающей «укрупненному» разбиению, относительно которого решается задача (8).
4. Определение обменной матрицы и поля скоростей среды. Метод последовательного усреднения
Нахождение обменной матрицы. Предположим, что область D, имеющая форму параллелепипеда в Rm, разбита на ячейки, которые мы для простоты считаем одинаковыми множествами: отрезками при m=1, квадратами при m=2 и кубами при m=3. В дальнейшем считается, что ячейки возникают в результате разрезания области точками (m=1), двумя семействами параллельных сторонам D прямых (m=2) или тремя аналогичными семействами плоскостей (m=3). Ребро ячейки обозначим через h. Эти ограничения на форму ячеек не обязательны, но сильно упрощают итоговые формулы.
Искомые матричные элементы aij=(mesDi)-1(Gej,ei). Их определение равносильно нахождению функций wj=Gej, которые являются решениями уравнений -cDw+(|Ñj|2/(4c))×w=ej при краевых условиях (6). Зафиксируем j и выпишем приближенные соотношения между элементами j-ого столбца матрицы aij. Каждый из этих элементов отвечает своей ячейке Di. Ячейки разделим на внутренние и граничные. Первые не примыкают к границе, а относительно вторых мы предполагаем, что, по крайней мере, одна из «граней» находится на границе области. Рассмотрим вначале внутреннюю ячейку, в которой находится источник с плотностью . К данной ячейке примыкают 2m соседних ячеек, имеющих с ней общие грани. Очевидно, что
. (11)
Здесь ¶kDj - граница с k-ой соседней ячейкой, а w является решением уравнения -cDw+(|Ñj|2/(4c))×w=ej при краевых условиях (6). Обозначим через - интегральное среднее по k-ой грани от
, матричный элемент для внутренней ячейки с источником через aив, а матричные элементы, отвечающие соседним ячейкам - через a(k). Получаем следующую формулу:
, при h ® 0.
Подставляя последнее выражение в (11) и обозначая значение ½Ñj½2 в центре выбранной ячейки через , получим равенство
Из последнего равенства получим
Обозначим матричный элемент, для внутренней ячейки, свободной от источника через acв. Аналогично предыдущему получим
Предположим, что ячейка является граничной. Соответствующий матричный элемент обозначим через aг. У такой ячейки в зависимости от m, количество s граней, примыкающих к границе, может быть равно 1, 2, или 3. Для граней, примыкающих к соседним ячейкам, как и раньше, считаем справедливым равенство
а для l-ой грани, примыкающей к границе, можно считать выполненным равенство
где sl - значение функции , определенной в (4), в центре l-ой грани.
Из (11), как и ранее получим для граничной клетки с источником выражение:
,
а для свободной граничной клетки - выражение:
.
Отбрасывая члены порядка o(h2), полученные соотношения между элементами j-го столбца обменной матрицы можно записать в векторно-матричной форме , где T - линейный оператор усреднения,
- вектор, определяемый источником. Если сходится матричный ряд
In + T + T2 + ¼ + TN + ¼,
то можно записать
.
Это равенство показывает, что вектор можно приближенно заменить вектором
, (12)
причём допускаемая при этом относительная погрешность (измеряемая в специально выбранной норме) не превышает величины , где
- операторная норма матрицы T. Для операторной нормы матрицы T, при следующем выборе нормы вектора
, справедливо неравенство
,
где . Таким образом, обратимость матрицы In-T и сходимость соответствующего матричного ряда гарантированы при условии v0 > 0. Разрешимость системы
и сходимость описанного выше алгоритма можно установить и без этого условия. Результаты вычислений показывают быстрое убывание с ростом N погрешности в приближенной формуле (12).
Метод, сведения задачи к рекуррентной формуле (12), назовем методом последовательного усреднения.
Нахождение потенциала поля скоростей среды также можно произвести методом последовательного усреднения. потенциал является решением следующей краевой задачи
;
где положительные функции считаются известными и удовлетворяющими условию
, которое означает, что приток среды в область D равняется величине стока. Эта краевая задача имеет множество решений, отличающихся постоянным слагаемым. Для выделения единственного решения будем считать выполненным еще одно условие
Пусть разбиение области D такое же, как и раньше. Среднее значение потенциала поля скоростей среды по ячейке Dj обозначим
Далее, как и выше, можно написать
Пусть - интегральное среднее по k-ой грани от
. Рассмотрим вначале внутреннюю ячейку Dj. Обозначим среднее потенциала для этой ячейки
, а средние для соседних ячеек - через
. Можно считать выполненным равенство
при h®0.
Таким образом,
Для граничной ячейки среднее потенциала обозначим . Обозначим количество граней, примыкающих к границе области, через
, а количества граней, примыкающих к частям
соответственно через
. Принимая во внимание (6), получим для ячейки, примыкающей к
Здесь - средние величины плотностей входящего и выходящего потоков среды по l-ой грани ячейки с номером j.
Для ячейки, примыкающей к , справедливо равенство
Отбрасывая погрешность аппроксимации o(h2), получим векторно-матричное уравнение
(13)
где T - линейный оператор усреднения, определяемый равенством
в котором предполагается, что для внутренней ячейки . Компоненты вектора
сопоставляются ячейкам и отличны от нуля только для ячеек, примыкающих к
причем для этих ячеек
Отметим, что матрица является вырожденной, поскольку для вектора
с одинаковыми координатами
Поэтому при произвольном векторе
система (13) не имеет единственного решения. Отметим также, что
. Если принять во внимание условие
, то мы должны разыскивать решение системы (13) в подпространстве векторов
удовлетворяющих условию
, а в этом подпространстве существует единственное решение. Хотя это подпространство не является инвариантным для оператора T, можно показать, что итерационный процесс
сходится к решению системы (13).
5. Численная реализация алгоритмов. Результаты вычислительных экспериментов.
Для численного решения рассмотренной задачи создан программный комплекс HeatCore, написанный на языках C#/C++. Имеется возможность решения как двумерных, так и трёхмерных задач с возможностью задания физических характеристик среды, параметров сеток и количества источников тепла. Краевые условия и функции задаются в виде скриптов на языке C#. Визуализация результатов расчёта: положение и мощность источников, поле скоростей, температурное поле осуществляется с помощью графического модуля, использующего OpenGL. Задача линейного программирования решается классическим симплекс-методом с использованием искусственного базиса. HeatCore использует равномерные сетки и позволяет решать задачу оптимизации распределения источников тепла для простых областей: квадрат (m=2) и параллелепипед (m=3). Для решения дифференциальных уравнений, кроме приведённого метода последовательного усреднения, так же реализован набор конечно-объёмных схем и методов решения разрежённых СЛАУ.
На рис. 1 приведена блок-схема алгоритма решения m-мерной задачи для движущейся среды с использованием метода последовательного усреднения для нахождения обменной матрицы. На рис. 2 приведены графики стабилизации значений (Jn)min с ростом n, полученные в результате численных экспериментов в трёхмерном случае, соответственно для неподвижной и движущейся сред с различными значениями коэффициента температуропроводности χ.
|
Рис. 1. Блок-схема алгоритма решения m-мерной задачи
|
а) неподвижная среда б) движущаяся среда
Рис. 2. Зависимость значения (Jn)min от n
На рис. 3 приведён результат численного эксперимента для 3-мерной области D размером 5×4×3 м с принудительным источником холода. В качестве вещества, заполняющего область, был выбран воздух с коэффициентом теплопроводности ϰ = 0,0244 Вт/(м∙К), плотностью = 1,293 кг/м3 и удельной теплоёмкостью c = 1005 Дж/(кг∙К) (расчёт произведён в системе СИ). Границы области выполнены из кирпича с коэффициентом теплопроводности 0,5 Вт/(м∙К) толщиной 30 см. В правую грань области входит вещество со скоростью s(x,y,z) = 10-4 м/с (
имеет площадь 2 м2). Симметрично с левой стороны имеется участок стока
такой же площади с такой же скоростью вещества s(x,y,z) = -10-4 м/с. Минимальный и максимальный профили температур заданы функциями m(x,y,z) = 290 К, M(x,y,z) = 320 К, температура внешней среды T0 = 260 К. Область контроля занимает всю область за исключением границы
.
Рис. 3. Оптимальное расположение плотности источников тепла внутри параллелепипеда.
Разноцветными объёмами показана искомая оптимальная плотность источников тепла (красные) и сток (синий). Менее прозрачные объёмы соответствуют более мощным источникам. Источники распределяются вблизи места входа вещества , границы, а также около стока, таким образом его нейтрализуя. Суммарная интенсивность источников в системе СИ имеет значение (Jn)min = 432,739 Вт.
Для проверки эффективности был проведён следующий эксперимент. Было взято некоторое случайное распределение источников (рис. 4а) и решена прямая задача теплопроводности для той же самой области (рис. 3), но без стока тепла. Параметры среды были взяты следующие: = 2,216×10-5 м2/с,
= 1,7×10−6 м/с. На дальней плоскости имеется окно с большим значением
= 3,4×10−5 м/с. В результате было получено, что случайное распределение источников нагревает область контроля в диапазоне температур 292,8-364,8 К. Затем была решена прямая задача с m(x,y,z) = 292,8 К и M(x,y,z) = 364,8 К при
. В итоге получилось, что случайное распределение (рис. 4б) имеет суммарную интенсивность 0,02687 К∙м3/с, а оптимальное – 0,01476 К∙м3/с. Сэкономленная мощность при этом составляет около 45 %.
а) неоптимальное б) оптимальное
Рис. 4. Расположение источников в параллелепипеде
6. Заключение. Для достаточно точного вычисления обменной матрицы требуется численно решать семейство краевых задач с использованием тонких сеток. Авторы так и поступали в работах [4]-[6]. Однако с ростом n в конечномерной аппроксимации Z0(n) количество этих краевых задач, равное количеству ячеек Di, становиться большим, особенно в трехмерном случае. Это может приводить к невозможности завершить вычисления за разумное время. Поскольку величина n заранее не известна, возникает потребность в менее затратном способе нахождения обменной матрицы, который позволил бы производить надежные численные эксперименты. Метод последовательного усреднения и является таким способом. Численные эксперименты показывают, что при относительно большом значении отношения последовательность (Jn)min с ростом n сравнительно быстро стабилизируется вблизи g и можно обойтись небольшим значением n. С уменьшением величины
скорость стабилизации этой последовательности снижается. Поэтому нахождение квазирешения в широком диапазоне значений параметров модели может потребовать использования аппроксимации Z0(n) с большим значением n. При этом далее приходится решать задачу линейного программирования с n переменными, что тоже может привести к значительным затратам машинного времени.
Отметим, что при одних и тех же параметрах модели результаты численных экспериментов, произведенных с помощью решения семейства краевых задач и с использованием метода последовательного усреднения практически совпадают.
1. Лионс Ж.-Л. Оптимальное управление системами, описываемыми уравнениями с частными производными. М.: Мир, 1972, 412 с.
2. Федоренко Р.П. Приближенное решение задач оптимального управления. М.: Наука, 1978, 497 с.
3. Алифанов О.М. Обратные задачи теплообмена. М.: Машиностроение, 1988, 280 с.
4. Брусенцев А.Г., Осипов О.В. Приближенное решение задачи об оптимальном выборе источников тепла // Научные ведомости БелГУ. Математика. Физика. №5 (124). Выпуск 26. Белгород, 2012. С.60-69.
5. Осипов О.В. Оптимальное расположение источников тепла в неоднородной среде // Вестник Белгородского государственного технологического университета имени В.Г. Шухова. №1. 2013. С.154-158.
6. Брусенцев А.Г., Осипов О.В. Оптимальный выбор источников тепла при наличии конвекции // Научные ведомости БелГУ. Математика. Физика. №26 (169). Выпуск 33. Белгород, 2013. С.64-82.
7. Михлин С.Г. Вариационные методы в математической физике. М.: Наука, 1970.
