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

64 НАУЧНЫЕ ВЕДОМОСТИ ЩШ Серия: Математика. Физика. 2013. №26(169). Вып. 33 MS С 26АЗЗ

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

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

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

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

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

А.Г. Брусенцев, О.В. Осипов

Белгородский государственный технологический университет им. В.Г. Шухова, ул. Костюкова, 46, Белгород, 308012, Россия, e-mail: brusentsev@mail.ru, ov.osipov@gmail.com

Введём функцию и(Х) = Т(X)-Т0, где Т(X) — установившаяся температура в точке X ограниченной связной области Д а Т0 — температура окружающей среды. Эта функция должна удовлетворять уравнению

где % — коэффициент температуропроводности среды, //(X) — поле скоростей среды, которое предполагается подчиненным условию йгьи = 0, а f (X) — плотность источников тепла в области Д Для формулировки краевых условий разобьём границу области Д на три части дД = Г+ и Г0 и Г_, оде Г+ — часть границы, которая является входом среды в область Д Г_ — часть границы, являющейся выходом (стоком) среды, а Г0 -часть непроницаемой границы. Справедливы следующие соотношения:

где п — единичный вектор внешней нормали к границе дД Последнее из этих трех условий означает, что в область Д поступает вещество внешней среды с температурой Т0

через непроницаемую для среды границу, равный а(х) ■ и(Х) (X Е Г0), оде а(х) >

0 — коэффициент теплопередачи через Г0, Следующие краевые условия выражают тепловой баланс области Д с окружающей средой.

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

Однако вопрос о существовании и единственности решения краевой задачи (1), (2) сложен дня анализа. В частном случае, когда ноле скоростей среды является потенциальным, эту задачу можно преобразовать к более простому і;илу.

Обозначим через ^(Х) — потенциал поля скоростей среды, т.е. /?(Х) = У^(Х), При этом У/?(Х) = Д^(Х) = 0, (п, У<р) |г = 0. Введём новую неизвестную функцию эд(х) с

помощью равенства и = ^е^/(2х), Подставляя это выражение в (1) получим следующее уравнение

1. Формулировка и преобразование краевой задачи установившегося теплообмена

хДи -У(и ■ и) + f = 0 , х Є Б

(1)

(п,г7) |Го=0, (п,г7) 1 Г- > 0 > (п,г7) |Г+ < 0

ХД^ + (|У^|2/(4х))^ = fe ^/(2х) .

Действительно, учитывая, что Д^(Х) = 0, получим

хД(^/(2х)) - У(тв^/(2х)У^) + / = х((А^)е^/(2х) + 2(Ут, Уе^/(2х)) + тДе^/(2х))-- (У т, У ^/(2х) - т(Уе^/(2х), У р) - те^/(2х) Д^ + / = = х((Д^)е^/(2х) + (1/х)(Ут, У^/(2х) + (1/(4х2)НУ^|2е^/(2х) + тД^/(2х))-- (Ум, У^/(2х) - (1/(2х))^|У^|2е^/(2х) - те^/(2х)Д^ + / = = е^/(2х)(хДм - (|У^|2/(4х))м) + /.

Отсюда получаем уравнение (3). Краевые условия примут следующий вид

(х(п, Ум) + ам) |Го = 0,

(х('Я,у«о + -('Я,У(/?)«о |г_=

21 |

(х('Я,у«о - -(я,У(/?)«о |г+=

0,

0.

Введём функцию, заданную на границе дБ образом

а(х)/х, а(х) = <) (п, У^)/2х;

-(п, У^)/2х,

Г+ и Г0 и Г- области следующим

х Є го; X Є Г— X Є Г+.

Тогда функция т должна удовлетворять краевому условию

(4)

дт

дп

ат

0.

(5)

дБ

Оператор Ьт = -хДт + (|У^|2/(4х))т, действующий в пространстве Ь2(Б) на достаточно гладкие функции т(Х), подчиненные краевым условиям (5), является самосопряженным и положительно определенным, а поэтому имеет ограниченный обратный оператор, определенный на Ь2(Д). Отсюда следует, что краевая задача (3), (5) имеет

Ь

Положительная определенность этого оператора вытекает из неравенства Фридрихса (см., например, |4|, с. 129)

J |Уи|2дУ + J и2дз ^ § J и2дУ ,

Б дО Б

где и(X) — произвольная функцня из С\(В) П С(5), а § > 0 — независящая от и(Х) константа. Действительно, для функций и(Х) € С2(5) П С(5) справедливо равенство

{1ли,и0= / |Угу|2йУ- /«’|^ + /*(1^|7(2х))ЛУ.

Б

дБ

Б

Если и (X) удовлетворяет краевым условиям (5), то можно с читать, что при X € ди

дИ = —а(х)т2. где а(х) ^ оо = сош! > 0. Таким образом, справедливо неравен-

дп

ство

(Ьи,и) ^ J |Уи|2дУ + а0 J и)2дв .

Б дБ

из которого с учетом неравенства Фридрихса получим

(Ьт,т) ^ §шт(а0,1)||и||2 ,

что и означает положительную определенность оператора Ь,

Наконец, отметим, что потенциал поля скоростей тоже является решением краевой задачи, которая состоит из уравнения

= 0

и краевых условий

(п, |Го = 0 ,

(п, У^) |Г_ = (X), (6)

(п, У^ |Г+ = 52(Х) .

где положительные функции 81(Х),52(Х) считаются известными и удовлетворяющими условию

J в1(Х)дБ = J з2(Х)д$, (7)

дГ_ дГ+

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

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

^(Х)дУт = 0. (8)

Б

Наконец отметим, в основной краевой задаче (3), (5) функция а(Х) в (5) может быть записана в виде

'а(Х)/х, Х € Го ;

1

—31 (;г), х Е Г_ ;

а(Х)

2х

—з2(х), геГ+.

2Х

Будем считать, что функция <^(х) известна. Температурный режим мы задаем в некоторой подобласти і С і, которую в дальнейшем называем областью контроля температуры. Нас интересует оптимизационная задача

где М(Х),т(Х) — соответственно максимальный и минимальный профили температуры в области контроля.

Учтем описанное выше преобразование краевой задачи. Получим следующую формулировку оптимизационной задачи. Определить минимум функционала (9) при условиях вы,пол,нения равенств (3), (5) и неравенств

Обозначим также через С оператор, обратный по отношению к оператору Ьт =

-хДт + (|У^|2/(4х))т. Тогда, принимая во внимание (10), (11), получаем следующую оптимизационную задачу

В оптимизационной задаче (12) Б(5) является классом функций, среди которых разыскивается решение #(Х). Этот класс желательно выбирать как можно шире, поскольку при этом можно рассчитывать на уменьшение экстремального значения 7{$}, Самым широким из возможных классов Б(5) можно было бы считать пространство суммируемых функций Ь1(5), поскольку на нем естественно определяется функционал 7{$}. Однако на этом пространстве не определен оператор С и задача (12) становится сформулированной не корректно. Корректной эта формулировка будет, если в качестве Б(5) взять пространство функций интегрируемых с квадратом Ь2(5). Ниже мы считаем, что Б(5) С Ь2(5). В особо оговоренных случаях считается, что Б(5) = Ь2(5).

2. Формулировка оптимизационной задачи

I = / (х)^У ^ тіп , (9)

Б

при условиях выполнения равенств (1), (2), а также неравенств

М(X) - Т0 ^ и(Х) ^ т(Х) - Т0, при X Є і; /(X) ^ 0, при X Є і , (10)

(М(X) - То)е-^/(2х) ^т^) ^ (ш^ - То)е-^/(2х) , при X Є і; /(X) ^ 0, при X Є і.

(и)

Введем обозначения

д^) = /^)е ^/(2х); б^) = (ш^) - Т0)е ^(ж)/(2х) 02(X) = (М(X) - Т0)е-^х)/(2х).

Б

б2^) ^ Сд^) ^ б^^), при X Є і; д^) Є 5(і),#^) ^ 0, при X Є і.

(12)

Если найдено точное или приближенное решение #0(Х) задачи (12)), то разыскиваемая оптимальная плотность источников те ила равна

/ш1п(Х) = в^(ж)/(2х)^о(Х). (13)

Замечание. Возможны и другие формулировки оптимизационной задачи. Иногда ограничения сверху в (12) отсутствуют или не существенны. Такую оптимизационную задачу назовем односторонней. В задаче (12) возможно появление дополнительных условий, которые могут привести к некоторым новым её модификациям. Одна из естественных модификаций состоит в требовании невозможности расположения источников тепла в некоторой части области 5, т.е. возникает дополнительное требование: #(Х) = 0 при Х € 50 С 5. Такую модификацию назовем задачей с ограничениями на, локализацию источников. Еще одна модификация связана с присутствием некоторых фиксированных источников до оптимизации. При этом функция д(Х) в (12) состоит из двух слагаемых, одно из которых известная функция, а второе подлежит определению. Такую задачу назовем задачей с фтксированными источниками. По форме она мало отличается от задачи (12). К функционалу 7{$} добавляется постоянное слагаемое, а функции 01(Х),02(Х) изменяются на однозначно определяемые слагаемые. Возможны различные комбинации рассмотренных выше модификаций задачи (12).

Мы не обсуждаем здесь вопрос о существовании точного решения задачи (12) и её модификаций. Вообще говоря, подобные задачи не обладают свойством единственности решения (см. |3|, замечание 1) Ответы на эти вопросы сильно зависят от выбора класса функций Б(5). Для практических нужд во многих случаях достаточно научиться находить так называемое квазирешение. Введем обозначение

7 = Б 7 Ы, (14)

деЯ(Б)

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

Определение 1. Квазирешеипем оптимизационной задачи (12) или её модификаций относительно класса Б(5) при данном допуске £ > 0 назовем такую функцию ^0(Х) € Б(5), удовлетворяющую ограничениям задачи, для которой выполняется неравенство

7Ы ^ 7 + £-

£

Б(5).

3. Конечномерная аппроксимация оптимизационной задачи

Дня приближенного нахождения квазирешения построим конечномерную аппроксимацию задачи (12) в виде задачи линейного программирования. Разобьём область 5 на

п частей , Определим подпространство Бп(Б) С Б(Б) кусочно-постоянных

функций вида д(Х) = ду, X € Б у (^ = 1,..., п). Введем в Бп(Б) базис, состоящий из функ-

ций бу (X) = 1, X € Бу, и бу (X) = 0 X $ Б у. Тощ а д(Х) = X] дубу (X). Разбиение обл асти Б мы

3 = 1

считаем и разбиением области Б С Б, т-е- ПРИ некотором натуральном р справедливо

а*у = (шее Б;) 1(Сбу, бг), а = (шее Бг) 1(01,б»), 6* = (швв Б*) 1(^2,бг),

где (•, •) — скалярное произведение в Ь2(Б), Заменяя класс функций Б (Б) подпространством Бп(Б), умножая ск^тарно ограничения на базисные функции бг(X), получаем конечномерную аппроксимацию задачи (12)

Здесь су = / б^(ж)/(2х)^У, Самым трудным при построении задачи (15) является на-

хождение матрицы агу = (шее Бг)-1(Обу, вг), поскольку оператор О явно не задан. В дальнейшем эта матрица называется обменной ,матрицей, Её определение равносильно нахождению функций Шу = Обу которые являются решениями уравнений —хДш + (|У^|2/(4х))ш = бу при краевых условиях (5). Вопрос о нахождении обменной матрицы мы рассмотрим ниже.

Так же как и в работе [3], задачу (15) мы обозначим через ^0(п) и будем называть конечномерной аппроксимацией задачи (12). Если задача (12) является односторонней, то и в задаче (15) тоже будут отсутствовать ограничения сверху Ьг = то. Такую задачу будем обозначать ^+(п). Аналогично можно получить конечномерные аппроксимации других модификаций задачи (12). Так дня задачи с ограничениями па локализацию источников в конечномерной аппроксимации уменьшается количество переменных. Исчезают некоторые слагаемые в целевой функции и ограничениях. Количество ограпи-

р

зации получим задачу вида (15) с другими значениями чисел аг, Ьг. В целевой функции появится свободный член.

п

приближенное квазирешепие задачи (12) или её модификаций.

4. Регулярность конечномерной аппроксимации по функционалу

Минимальное значение целевой функции задачи (15) или конечномерной аппроксимации любой другой модификации задачи (12) будем обозначать через (/п)т;п.

П

Р

равенство Б = и Б Б(Б) Од

г=1

п

п

(15)

Определение 2. Конечномерную аппроксимацию назовем регулярной но функционалу в классе Б (Б), если справедливо неравенство

где 7 определено в (14),

При наличии регулярности решение конечномерной аппроксимации задачи при достаточно малом шах(&аш Б) можно считать приближенным квазирешением исходной

г

задачи. Действительно, после дискретизации система ограничений означает, что неравенства в исходной задаче практически выполнены, т. к. они удовлетворяются в среднем по Б. При этом значение (3п)тт при достаточно больших п не превосходит 7 + £, Ниже мы докажем регулярность конечномерных аппроксимаций ^+(п) и Дз(п) по функционалу в классе Б (Б) = £2(Б).

В дальнейшем для аппроксимации ^о(п) будем считать выполненным следующее Условие А. Начиная с некоторого номера у задачи ^0(п) существует хотя бы одна точка минимума, для которой правые неравенства в (15) являются строгими.

Мы не станем обсуждать вопрос о том, когда это условие выполнено. Отметим лишь, что аппроксимацию ^+(п) можно всегда рассматривать как задачу ^0(п) (Ь = то), удовлетворяющую условию А.

Установим несколько вспомогательных утверждений. Вместе с последовательностью задач ^0(п) мы будем рассматривать задачи Z(^ (п), построенные заменой а и Ь на а — 6 иЬг + 6 при некотором 6 > 0,

Лемма 1. Если задача Z0(n) имеет решение, то при любом 6 > 0 справедливы неравенства

где у* (г= 1, 2/г) — координаты какой-нибудь точки максимума задачи, двойственной но отношению к Z0(n), а /п(6) — минимальное значение целевой функции задачи Z<^ (п).

□ Поскольку точка минимум а задачи Z0(n) удовлетворяет системе ограничений задачи Z(^(п), то справедливо неравенство (3п)тт ^ 1п(6), Целевая функция задачи, двойственной по отношению к задаче Z(^(п), имеет вид

Согласно первой теореме двойственности, если обозначить через 3^(у) целевую функцию задачи, двойственной к Z0(n), а терез у0 — точку максимума задачи двойственной к Z(^ (п), получим

2п

i=1

п

п

2п

i=l

i=l

i=1

2п

2п

2п

1.(6) = (у0) — 6 £ у0 > (У •) — 6 £ у* = (3п)т,п — 6 £ у* ,

i=1

i=1

i=1

что завершает доказательство. ■

Лемма 2. Пусть (3п)тт и 1п($) — минимальные значения целевых функций задач ^о(п) и Z(^ (п) соответственно. Если выполнено условие А, то справедливы неравенства

(3п)тт ^ 1га(^) ^ (3п)тт(1 ^ ' (ш-1Н аг) ) •

г

□ В силу условия А для точки минимума задачи Z0(n) ограничения сверху в (14) являются строгими неравенствами, а это, в силу условий дополнительной нежёсткости, означает, что в точке максимума двойственной задачи у* = 0 (г = п, 2/г). Поэтому

п 2п

Штт = (У Ч = X! агУг* ^ (™П аг) X! Уг* •

г=1 г=1

Отсюда и из леммы 1 вытекает справедливость утверждения. ■

В следующей лемме фигурирует линейный функционал вида

3Р{д} = J р(Х)д(Х)^,

Б

определенный на функциях д(Х) € £2(Б). Весовая функция р считается удовлетворяющей неравенствам С1 ^ р(Х) ^ С2 (X € Б) с положительными константами С'1,С'2.

Лемма 3. Для произвольного числа а > 0 и любой функции д(Х) € Ь2(Б) найдется такое число £ > 0, что при выполнении для разбиения области Б неравенства

П

шах(&ашБз-) < £ можно выбрать кусочно-постоянную функцию д(Х) = ^ дз-вз-(X), для

у=1

которой выполнены условия 3р{д} = 3р{д}, ||д — д|| < а гДе II ■ II _ норма в Ь2(Б). При этом если почти всюду в Б д(Х) ^ 0, го дз- тоже можно выбрать неотрицательными.

□ Для любой функции д(Х) € Ь2 (Б) можно найти такую функцию ^(Х) € С0(С), что ||^ — д|| < аС^/^СУ2), Ввиду равномерной неирерывности ^(Х) существует такое

£ > 0, что для разбиения области Б при условии шах(^гатБз-) <£ можно найти кусочно-

з

постоянную функцию д^Х), для которой |д1 — < аС1/2/(2С11/2) или |д1 — д|| <

аС^/С^2, Зафиксируем это разбиение области Б и рассмотрим задачу на минимум для функции п переменных

П р

ф(д) = 11р1/2(д — д)|2 = ^(дг2 ■ ( / Р(Х)^) — 2дг ■ (рд,вг)) +

г=1 Б

В стационарной точке справедливы равенства д» • (/ р(х)с1У) — (рд, е») = 0, г = 1, /г. Вы-

Бг

бирая дг = (рд,вг)// р(Х)а!У, получим кусочно-постоянную функцию д(Х), для которой 3{д} = 3{д}. С другой стороны

С,!1(д — д)||2 < ф(д) « Ф(д1) = ||р1/2(д! - д)||2 < С2|д, — д||2 < С,а2 •

НАУЧНЫЕ ВЕДОМОСТИ К,Д Серия: Математика. Физика. 2013. №26(169). Вып. 33 73 Поэтому справедливо неравенство ||д — д|| < а. Ш

Лемма 4. Если т < 4 и выбрано разбиение области И такое, что Б = У Б* С Б,

то для любого элемента разбиения Б* С Б, і = 1 ^ р и любой функции д(Х) Є Ь2(Б) выполнено неравенство

где константа С> 0 не зависит от I, выбора разбиения области Б и функции д.

(2)

□ Отметим, что функция Сд принадлежит Соболевскому пространству Ж (Б) и является непрерывной. Поэтому по теореме о среднем интегрального исчисления получаем

Если учесть еще одно неравенство Фридрихса |6, с. 196, предложение 3.6.111, а также положительную определенность оператора Ь получим

Теорема. Пусть т < 4 и в задаче (12) Б С Б, а также при всех Х € Б и некотором ^0 > 0 выполнено неравенство 01(Х) ^ ^0. Пусть также последовательность разбиений области Б удовлетворяет условию Иш шах(^гатБз-) = 0. Тогда, если выполнено усло-

вие А, то соответствующая последовательность конечномерных задач ^з(п) регулярна по функционалу в классе Ь2(Б).

□ Пусть а — произвольное положительное число. Выберем функцию 0 ^ д(Х) Є Ь2(Б), удовлетворяющую условиям (12), и такую, что 7{д} — 7 < а, оде 7 — число, определенное в (14). Согласно лемме 3, вследствие условий теоремы, дня всякого разбиения области Б, связанного с задачей ^з(п) при достаточно большом и, найдется кусочно-постоянная функция д(Х), для которой выполнены условия 7{д} = 7{д}, ||д — д|| < а. Эта функция может не удовлетворять неравенствам в (12). Однако её координаты д* удовлетворяют системе ограничений задачи Z(^ (и) при некотором 6 > 0 и достаточно большом и. Действительно, при Х Є Б справедливы неравенства 02(Х) ^ Сд = Сд + С(д — д) ^ ^(Х), Умножим скалярно в Ь2 (Б) эти неравенств а на е*(Х) и разделим на шее Б* при і = 1, 2,... ,р. Вследствие леммы 4 получим

р

і=1

|(шевБ*) 1(Сд,е*)| ^ С||д|и2(д)

По теореме вложения |5, с. 74, теорема 1| справедливо неравенство

тах |Сд| ^ Ц^дЦ^)^

п

(шее б і) 1(сд,ег) = ^ а, д, ^ щ — (шее б і) 1|(с(д — д), е*) | ^

і = 1

(шея Д) —з' < Ъг + (шея Дг) 1|(С(^ - £),е*)| <

3=1

< Ъг + С||д - £|| < Ъг + Са .

Тем самым показано, что д(Х) удовлетворяет системе ограничений задачи ^(п) с 5 = Са. Отсюда, если учесть лемму 2, получим

JШ ^ 1п(5) ^ (Д)т1п(1 - Ск- 1а).

Таким образом, для любого сколь угодно малого а> 0 найдется такой номер Ма, что при п ^ Ма выполнено неравенство (Д)т1п ^ (т+а)(1 - Ск- 1а)- 1, которое и доказывает теорему. ■

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

5. Разностная схема для вычисления поля скоростей

Построим иптегро-иптерполяциоппым методом |7| разностную схему для расчёта поля скоростей в двумерной и трёхмерной областях. В двумерном случае в качестве области Д возьмём прямоугольник, в 3-мерном — параллелепипед. Сетку, на которой будем искать решение, возьмём равномерную. В двумерном случае сетка имеет вид (рис. 1). Обязательно должны быть узлы в углах области и па гранях.

В т-мерном случае для области Д = (аа ^ ха ^ Ъа, а = 1,... ,т) сетка имеет вид:

{хг (а1 + ^1 ^15 ... 5 ат + imhm) , ^ (^15 ... 5 ^т) 5 ° . . . 5 ha (Ъа аа) /^а} .

Окружим каждый внутренний узел хг (га = 15... 5 N — 1) элементарным объёмом Ег = {аа + 1аЬа - Д/2 < Ха < аа + ІаД + Д/2} .

Для приграничных узлов хг(га = °5 iа = Д, - 1) элементарные объёмы ограничиваются границей области Д как показано та рис. 1 (объёмы До,о), Д0,3))-

Рис. 1. Разбиение области Д сеткой размером 4 х 4.

Обозначим также , £“+ - поверхности соприкосновения объёма Е" с соседними по пространственному измерению ха объёмами Е|а+, г а+ = (г1;..., га + 1,..., гт) и Е^а-, га- = (г1,...,га — 1,..., гт) или с границей области Е, если Е" — приграничный объём и он какими-то из граней касается границы. Р“-, Ра+ — потоки через соответствующие поверхности. Преобразуя объёмный интеграл по формуле Гауееа-Остроградского в поверхностный, получаем:

/т т я а я

=£(рг+рп = £ £™+

а=1 а=1 ±_ „а

( \

дф

-^-сІБ =0. 16)

дп

а— па +

\?г *?+ у

Аппроксимируем функцию (внутри области Е" (га = 1,..., N — 1) сеточной функцией ( и запишем предыдущее уравнение в виде:

т

ка * к

Учитывая, что шее £“+ = шее £“ = шее Е|/ка, получим

J А^Ут = ^ ----^2г + ^га+^ шее Е? + 0(Л2) = О . (17)

Ег

Если какая-либо поверхность объёма Е| касается границы (га = 0, га = У«), то в уравнении (17) соответствующее значение потока формируется исходя из условия (6). Например, если вещество входит в область Е через границу Ж"-, то:

РГ= /^»= [»№>■

Уравнение (8) аппроксимируется уравнением

п N Мт

/ ^ ^ ((/3*шее Е*) + 0(к2) = 0. (18)

^ *1=0 *т=0

Составленные для каждого объёма Е* уравнения (17) образуют систему линейных уравнений с неизвестными ( общим количеством (Ж1 + 1) ■.. .■ (Ат +1), имеющую решение с точностью до произвольной постоянной. Если к ней добавить уравнение (18), то система будет иметь единственное решение, при условии, если выполняется условие (7). Результатом её решения будет сеточная функция потенциала ф.

а—

S

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

Если область Е какой-то поверхностью касается границы (например, по направлению против оси ха), то, согласно условию (5), поток через неё равен:

Соответствующие слагаемые в этом случае также нужно заменить в уравнении (19). Составленные для каждого объёма Е уравнения (19) образуют самосопряжённую систему линейных уравнений с неизвестными гй|, имеющую единственное решение. Система имеет п правых частей. Её решение даёт набор сеточных функций г- (X) (^ = 1,..., п), каждая из которых является приближённым решением (3), (5) с f (X) = е- (X), и можно считать, что

Далее, для вычисления элементов а- используется численное интегрирование методом трапеций. Подобласти Д- выберем следующим образом: разрежем область Б по каждому измерению па равные части. Таким образом, получим одинаковые но объёму подобласти, которые образуют также равномерную сетку ш/, Сетка ш должна быть чаще ш/ по каждому измерению в целое число раз п чтобы узлы сетки ш были в точках разрыва функций е— Если расчётная сетка ш крупная, то можно просто считать

6. Нахождение элементов обменной матрицы а,

+ ((|У(|2/(4Х)№ - у /е-^(2х))^т + 0(к2) = 0 . (19)

Се,

т

а, = т,(х)е^(ж?)/(2х)шее А, иначе а, = / т,(Х)е^(х)/(2х)^Ут.

ХєПі

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

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

Для численного решения задачи (15) при т =2, 3 разработан программный комплекс Не^Соге, написанный па языке Сф и не использующий дополнительных математических библиотек. Данный комплекс позволяет принимать все необходимые входные данные: константы, коэффициенты, параметры сеток, в том числе краевые условия, функции температур в виде скриптов па языке Сф. Для представления результатов разработаны графические модули дня вывода:

1. Графиков функций 2 переменных в виде поверхностей на 3-мерпой сцене (рис. 5);

2. Кусочно-постоянной функции плотности источников па 3-мерпой сцепе в виде объёмов различной прозрачности и цвета (рис. 3);

3. Тепловых карт и линий тока вещества (рис. 6).

На рис. 2 приведена общая блок-схема алгоритма решения задачи (15).

Теоретически в работе установлено, что при стремлении шах(&аш Д-) к нулю точное

-

значение (/га)тт стабилизируется около нижней грани целевого функцианала. Поэтому п мы выбираем эмпирически по признаку стабилизации (/га)т;п.

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

Эксперимент 1. Области Д, заполненной воздухом с коэффициентом температуропроводности х = 2, 216 х 10-5 м2/с зададим раз мер 5 х 5 х 3 м. Границу будем считать выполненной из материала с хк = 5, 2 х 10-7 м2/с толщиной 30 см (кирпичная стена), местами из материала Хс = 3, 4 х 10-7 м2/с толщиной 1 см (стекло). Краевые условия дня расчёта поля скоростей зададим следующим образом: в местах притока вещества Г+ з2(ж) = 0, 0002 м/с, в местах стока Г- ^(ж) = 0, 0002 м/с; размер расчётной сетки 15 х 15 х 15

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

(20)

5

Б

7. Результаты численных экспериментов

т(ж,у,г)

285 - (у - 4, 5) ■ 12, 5 К, у ^ 4, 5 и г ^ 1;

285 К, у ^ 4, 5 ил и г ^

и М(ж, у, г) = 310 , Функция т(ж, у, г) выбирается меньшей в месте входа вещества, потому что нагреть вещество в этом месте можно только слишком сильным источником, который также сильно разогреет всю область. Температуру внешней среды возьмём Т0 = 260 К. Далее будем считать область, где у ^ 4, 5 или г ^ 1 областью контроля. На рис. 3 изображён результат решения задачи: полупрозрачными объёмами показаны источники, прозрачность которых тем меньше и цвет ярче, чем больше интенсивность источника. Для данного расчёта результирующее значение (/га)т;п = 0, 018497 К ■ м3/с. Учитывая плотность воздуха р = 1, 293 кг/м3 и его удельную теплоёмкость св = 1005 Дж/(кг-К), получим мощность источников в более привычной физической системе единиц, равную 0, 018497 ■ 1, 293 ■ 1005 Вт и 24, 036 Вт.

Рис. 2. Блок-схема алгоритма решения ш-мерной задачи.

Замечание. Для данного эксперимента минимальное значение потенциала ^min ~

—0, 000199 м2/с наблюдается в точке входа вещества. При этом в правой части уравнения (3) / умножается на exp(1, 99 х 10-4/(2 • 2, 216 х 10-5)) ~ e4,49. Дальнейшее увеличение скорости приводит к появлению чисел больших порядков и ошибкам округления при использовании тина double. Например, если ускорить движение вещества в 10 раз, числа будут размера e44,9 ~ 3,16 х 1019, Стандартный вещественный тип данных double, как известно, держит только 15-16 значащих цифр, поэтому для подобных алгоритмов нужно использовать тины повышенной точности.

Рис. 3. Оптимальное расположение источников тепла / в параллелепипеде.

Эксперимент 2. При тех же входных данных и параметрах среды проведём вычисления с различным % и числом разбиений п. Частоту сетки п достаточно брать равной

2, поскольку увеличение п не оказывает существенного влияния на конечный результат. Увеличение % приводит к более быстрой по п стабилизации вели чины (/га)т;п (табл. 1).

Таблица 1

Зависимость (.1п)ттОТ х И П

п \ х 2 х 11) :'.\г с 4x10 5м2/е 8 х 11) :'.\г с 1, 6 х 10_4м2/е 3, 2 х 11) '.\г с

4x4x4 0,02175 0,01978 0,01723 0,01514 0,01379

5 х 5 х 5 0,01872 0,01641 0,01463 0,01349

6x6x6 0,01918 0,01781 0,01587 0,01430 0,01331

7x7x7 0,01888 0,01747 0,01571 0,01426 0,01330

8x8x8 0,01880 0,01758 0,01585 0,01428 0,01330

9x9x9 0,01893 0,01759 0,01575 0,01421 0,01328

10x10x10 0,01863 0,01734 0,01556 0,01413 0,01322

11x11x11 0,01817

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

Таблица 2

Зависимость (7п)тт ОТ ^(Х) И П

50 п 10_5м/с ю X О 1 СЛ £ о' X О 1 СЛ £ о' 8х10“5 м/с 1,6х10“4 м/с 3,2 х 10-4 м /с

4x4x4 0,00521 0,00581 0,00723 0,01058 0,01790 0,03141

5x5x5 0,00516 0,00570 0,00696 0,01003 0,01694 —

бхбхб 0,00515 0,00568 0,00689 0,00975 0,01590 0,02791

7x7x7 0,00511 0,00562 0,00680 0,00956 0,01564 0,02778

8x8x8 0,00509 0,00559 0,00676 0,00956 0,01562 0,02728

9x9x9 0,00511 0,00561 0,00678 0,00957 0,01567 —

10x10x10 0,00509 0,00558 0,00672 0,00944 0,01547 0,02734

Эксперимент 4. Для проверки эффективности использования разработанного алгоритма проведём следующий эксперимент. При тех же входных данных и параметрах среды расположим источники случайным образом (Рис. 4а). Данные источники суммарной мощностью 0, 094491 К-м3/с обеспечивают температурный коридор в области контроля в диапазоне 278-578,4 К. Эти значения возьмём в качестве входных данных для решения задачи оптимизации (15). Пусть в области контроля т(х,у,г) = 278 К, М(х, у, г) = 578,4 К. Оптимальное распределение источников (Рис. 46) при этом даёт суммарную мощность 0, 01335 К-м3/с, что экономит больше 80% тепловой энергии.

а) неоптимальное б> оптимальное

Рис. 4. Расположение источников в параллелепипеде.

Эксперимент 5. Для т =2 результат решения задачи (плотность источников) представлен в виде столбиков, и результирующая температура в виде поверхности (рис. 5). На сторонах расположены места стоков - белые параллелограммы. Приток находится па дальней стороне около самого сильного источника. Параметры данного эксперимента приведены на рис. 5.

. г, К

3O3

у, м 3OOi

Г (сток) Si(x,y) = 1

(■JL„ = 257,744 К-м2/с

lOxlO TO = 26O К X = 5 м2/с П = 4

m(x,y) = 3OO K М(х,у) = 3O5 K

a(x, O) = l м/с

f

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

Г+ (приток) s2(x,y) = І м/с

Рис. 6. Линии тока вещества v.

Литература

1. Брусенцов А.Г., Брусенцова B.C. Задача об оптимальном выборе источников тепла /7 Сб. трудов XXIII международной конференции «Математические методы в технике и технологиях». 2010. 2. С.43-46.

2. Брусенцов А.Г., Осипов О.В. Численное исследование задачи об оптимальном выборе источников тепла /7 Сборник трудов XXIV международной конференции «Математические методы в технике и технологиях» (ММТТ-24), Саратов. 2011. 2. С.33-34.

3. Брусенцов А.Г., Осипов О.В. Приближенное решение задачи об оптимальном выборе

источников тепла /7 Научные ведомости БелГУ. Математика. Физика. 2012. №5

(124);Вып.26. С.60-69.

4. Михлин С.Г. Вариационные методы в математической физике / М.: Наука, 1970.

5. Соболев С.Л. Некоторые применения функциональнох'о анализа в математической физике / М.: Наука, 1968.

6. Нарасимхан Р. Анализ на действительных и комплексных мнох’ообразиях / Издательство «Платон», 1997.

7. Самарский А.А. Теория разностных схем / М.: Главная ред. физ.-мат. лит. изд-ва «Наука», 1977. 656 с.

OPTIMAL CHOICE OF HEAT SOURCES AT THE PRESENCE OF CONVECTION

A.G. Brusentsev, O.V. Osipov

Belgorod State Technological University named V.G. Shuchov,

Kostyukova St., 46, Belgorod 308012, Russia, e-mail: brusentsev@mail.ru, ov.osipov@gmail.com

Abstract. It is proposed and justified the new method of numerical solution of the optimal choice of heat source in moving medium. It consists in reducing the original problem to linear programming. In the course of numerical simulations using finite-diffcrencc schemes obtained numerical results in three-dimensional region with parameters closed to real ones at small velocity of fluid flow.

Key words: moving medium, optimal control, elliptic boundary value problems, finite-dimensional approximation, inverse problem of heat conduction, convection, density of heat, heat balance method, simplex method.