Журнал вычислительной математики и математической физики, 2021, T. 61, № 10, стр. 1594-1609

Определение коэффициента теплопроводности вещества по тепловому потоку на поверхности трехмерного тела

А. Ф. Албу¹, В. И. Зубов^1, *

¹ ФИЦ ИУ РАН
119333 Москва, ул. Вавилова, 40, Россия

Поступила в редакцию 30.01.2021
После доработки 30.01.2021
Принята к публикации 09.06.2021

DOI: 10.31857/S0044466921100033

Аннотация

Рассматривается и исследуется обратная задача определения зависящего от температуры коэффициента теплопроводности вещества по известному тепловому потоку на границе тела. Рассмотрение проводится на основе первой краевой задачи для трехмерного нестационарного уравнения теплопроводности в параллелепипеде. Обратная коэффициентная задача сводится к вариационной задаче и решается численно с помощью градиентных методов минимизации функционала. В качестве целевого функционала выбрано среднеквадратичное отклонение рассчитываемого теплового потока на поверхности тела от экспериментально полученного потока. На примере ряда нелинейных задач, коэффициенты которых зависят от температуры, показаны работоспособность и эффективность предложенного подхода. Библ. 18. Фиг. 10.

Ключевые слова: обратные коэффициентные задачи, нелинейные задачи, трехмерное уравнение теплопроводности, оптимальное управление, численные методы оптимизации, схемы переменных направлений.

ВВЕДЕНИЕ

Обратные задачи всегда были в поле внимания исследователей. Однако в последнее время наблюдается повышенный интерес к этим задачам. Причиной повышенного интереса является возросшая потребность в создании новых, современных материалов и изучении их свойств. Нередко некоторые характеристики вновь созданных материалов оказываются заранее неизвестными. Одним из способов их определения является решение обратных задач (см., например, [1], [2]).

Одной из важных и сложных обратных задач является задача определения зависимости коэффициента теплопроводности вещества от температуры по результатам экспериментального наблюдения за динамикой температурного поля в исследуемом веществе. Эта обратная коэффициентная задача исследовалась в работах [3]–[9]. В этих работах рассмотрение проводилось на основе первой краевой задачи для одномерного и двумерного нестационарного уравнения теплопроводности. В качестве исследуемых объектов там выступали стержень и прямоугольная пластина. Результаты проведенных исследований показали, что в двумерном случае получить решение обратной задачи существенно сложнее, чем в одномерном случае. Связано это не только с увеличивающимися затратами ресурсов компьютера, но и с особенностями распределения экспериментальных данных. Для получения решения обратной задачи в двумерном случае пришлось модернизировать численный алгоритм.

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

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

Для успешного численного решения обратной коэффициентной задачи в трехмерном случае необходимо использовать эффективные разностные схемы, аппроксимирующие прямую и сопряженную задачи. В работе [13] проведены исследования, касающиеся выбора разностной схемы для аппроксимации уравнения диффузии тепла при решении обратной коэффициентной задачи в трехмерной постановке. В [13] на примере ряда нелинейных задач для трехмерного уравнения теплопроводности, коэффициенты которого зависят от температуры, проводился сравнительный анализ нескольких схем переменных направлений: локально-одномерной схемы (см. [14]), схемы Дугласа–Рекфорда (см. [15]) и схемы Писмена–Рекфорда (см. [16]). При сравнении методов принимались во внимание точность получаемого решения и время достижения требуемой точности на компьютере. Результаты проведенных численных экспериментов показали, что локально-одномерная схема для таких задач является наиболее эффективной и наименее “капризной”.

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

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

1. ПОСТАНОВКА ЗАДАЧИ

Пусть образец исследуемого вещества имеет форму прямого параллелепипеда длины $R$, ширины $E$ и высоты G. Температура $T$ этого параллелепипеда в начальный момент времени известна. Известно также, как изменяется во времени температура параллелепипеда на его гранях. Для математического описания процесса теплопроводности в параллелепипеде введем декартовы координаты $x$, $y$ и $z$. Точки $s = (x,y,z)$ параллелепипеда образуют область $Q = \left\{ {(0,R) \times (0,E) \times (0,G)} \right\}$ с границей $\Gamma = \partial Q$. Распределение температурного поля в параллелепипеде в каждый момент времени описывается решением следующей начально-краевой (смешанной) задачи:

(1.1)

$C(s)\frac{{\partial T(s,t)}}{{\partial t}} = {{\operatorname{div} }_{s}}\left( {K(T(s,t)) \cdot {{\nabla }_{s}}T(s,t)} \right),\quad (s,t) \in \left\{ {Q \times (0,\Theta ]} \right\},$

(1.2)

$T(s,0) = {{w}_{0}}(s),\quad s \in \bar {Q},$

(1.3)

$T(s,t) = {{w}_{\Gamma }}(s,t),\quad s \in \Gamma ,\quad 0 \leqslant t \leqslant \Theta .$

Здесь $t$ – время; $T(s,t) \equiv T(x,y,z,t)$ – температура вещества в точке $s$ с координатами $(x,y,z)$ в момент времени $t$; $\;C(s)$ – объемная теплоемкость материала; $K(T)$ – коэффициент теплопроводности; ${{w}_{0}}(s)$ – заданная температура параллелепипеда в начальный момент времени $t = 0$; ${{w}_{\Gamma }}(s,t)$ – заданная температура на границе области. Объемная теплоемкость $\;C(s)$ считается известной функцией координат. При заданном коэффициенте теплопроводности $K(T)$ можно найти распределение температуры $T(s,t)$ в любой точке области $Q \times (0,\Theta )$, решив прямую задачу (1.1)–(1.3).

Обратная коэффициентная задача сводится к следующей вариационной задаче: требуется найти такую зависимость коэффициента теплопроводности вещества $K(T)$ от температуры $T$, при которой поток тепла $\left( { - K(T(s,t))\frac{{\partial T(s,t)}}{{\partial{ \bar {n}}}}} \right)$ на границе объекта, полученный в результате решения прямой задачи (1.1)–(1.3), мало отличается от потока тепла $P(s,t)$, полученного экспериментально. В качестве меры отклонения этих функций рассматривается величина

(1.4)

$\Phi (K(T)) = \int\limits_0^\Theta {\mathop{{\int\!\!\!\!\!\int}\mkern-21mu \bigcirc}\limits_\Gamma {\beta (s(\Gamma )){{{\left[ {\left( { - K(T(s(\Gamma ),t)){{{\left. {\frac{{\partial T(s,t)}}{{\partial{ \bar {n}}}}} \right|}}_{{s \in \Gamma }}}} \right) - P(s(\Gamma ),t)} \right]}}^{2}}d\Gamma } dt} .$

Здесь $\beta (s(\Gamma )) \geqslant 0$ – заданная весовая функция; $P(s(\Gamma ),t)$ – заданный тепловой поток на границе $\Gamma $ области $Q$; $\partial T{\text{/}}\partial{ \bar {n}}$ является производной от температуры по направлению внешней нормали к границе области.

Таким образом, обратная задача идентификации функции $K(T)$ сводится к следующей задаче оптимального управления: определить оптимальное управление $K(T)$ и соответствующее решение $T(s,t)$ задачи (1.1)–(1.3), при котором целевой функционал (1.4) достигает минимального значения.

2. ВЫЧИСЛЕНИЕ ГРАДИЕНТА ЦЕЛЕВОГО ФУНКЦИОНАЛА

Вариационная задача (1.1)–(1.4), к которой сводится обратная коэффициентная задача, решалась в данной работе градиентными методами минимизации функционала. В работе [17] было получено аналитическое выражение для градиента целевого функционала в общей, $n$-мерной постановке сформулированной задачи оптимального управления. Однако использовать аналитическое выражение для градиента функционала при численном решении задачи практически невозможно в силу его сложности. Поэтому в предлагаемом алгоритме, как и ранее, для вычисления градиента функционала использовалась методология быстрого автоматического дифференцирования.

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

Для дискретизации задачи в области $Q \times (0,\Theta )$ вводятся временнáя и пространственные сетки. Временнáя сетка задается набором узловых значений $\{ {{t}^{j}}\} _{{j = 0}}^{J}$, ${{t}^{0}} = 0$, ${{t}^{J}} = \Theta $ с равномерным шагом $\tau = \Theta {\text{/}}J$.

Вводятся также две пространственные сетки: основная и вспомогательная. Для построения основной пространственной сетки на отрезке $[0,R]$ выбирается система “опорных” точек $\{ {{x}_{n}}\} _{{n = 0}}^{N}$ так, что ${{x}_{0}} = 0$, ${{x}_{N}} = R$ и ${{x}_{n}} < {{x}_{{n + 1}}}$ для всех $0 \leqslant n < N$. При этом $h_{n}^{x}$ – расстояние между опорными точками ${{x}_{n}}$ и ${{x}_{{n + 1}}}$, т.е. $h_{n}^{x} = {{x}_{{n + 1}}} - {{x}_{n}}$, $n = \overline {0,\;N - 1} $. Аналогично на отрезках $[0,E]$ и $[0,G]$ выбирается система “опорных” точек $\{ {{y}_{i}}\} _{{i = 0}}^{I}$ и $\{ {{z}_{l}}\} _{{l = 0}}^{L}$ так, что ${{y}_{0}} = 0$, ${{y}_{I}} = L$, ${{y}_{i}} < {{y}_{{i + 1}}}$ для всех $0 \leqslant i < I$ и ${{z}_{0}} = 0$, ${{z}_{L}} = G$, ${{z}_{l}} < {{z}_{{l + 1}}}$ для всех $0 \leqslant l < L$ соответственно. При этом $h_{i}^{y} = {{y}_{{i + 1}}} - {{y}_{i}}$, $i = \overline {0,\;I - 1} $, $h_{l}^{z} = {{z}_{{l + 1}}} - {{z}_{l}}$, $l = \overline {0,\;L - 1} $. Точки основной сетки – совокупность точек $\{ {{x}_{n}},{{y}_{i}},{{z}_{l}}\} $, где $n = \overline {0,\;N} $, $i = \overline {0,\;I} $ и $l = \overline {0,\;L} $.

Вспомогательная сетка $\{ {{\tilde {x}}_{n}},{{\tilde {y}}_{i}},{{\tilde {z}}_{l}}\} $, $n = \overline {0,\;N + 1} $, $i = \overline {0,\;I + 1} $, $l = \overline {0,\;L + 1} $, строится похожим образом. Отличие состоит лишь в выборе опорных точек:

${{\tilde {x}}_{0}} = {{x}_{0}},\quad {{\tilde {x}}_{{N + 1}}} = {{x}_{N}},\quad {{\tilde {x}}_{n}} = {{x}_{{n - 1}}} + h_{{n - 1}}^{x}{\text{/}}2,\quad n = \overline {1,\;N} ,$

${{\tilde {y}}_{0}} = {{y}_{0}},\quad {{\tilde {y}}_{{I + 1}}} = {{y}_{I}},\quad {{\tilde {y}}_{i}} = {{y}_{{i - 1}}} + h_{{i - 1}}^{y}{\text{/}}2,\quad i = \overline {1,\;I} ,$

${{\tilde {z}}_{0}} = {{z}_{0}},\quad {{\tilde {z}}_{{L + 1}}} = {{z}_{L}},\quad {{\tilde {z}}_{l}} = {{z}_{{l - 1}}} + h_{{l - 1}}^{z}{\text{/}}2,\quad l = \overline {1,\;L} .$

Плоскости $x = {{\tilde {x}}_{n}}$, $y = {{\tilde {y}}_{i}}$, $z = {{\tilde {z}}_{l}}$ делят объект на параллелепипеды – так называемые расчетные ячейки. Расчетной ячейке будем приписывать индексы $(n,i,l)$, если вершина этой ячейки, наиближайшая к началу координат, совпадает с узловой точкой $({{\tilde {x}}_{n}},{{\tilde {y}}_{i}},{{\tilde {z}}_{l}})$. В каждом узле $({{x}_{n}},{{y}_{i}},{{z}_{l}},{{t}^{j}})$ расчетной области $Q \times (0,\Theta )$ все функции задаются своими значениями в этой точке (например, $T({{x}_{n}},{{y}_{i}},{{z}_{{l,}}}{{t}^{j}}) = T_{{nil}}^{j}$).

Обозначим через $\Gamma _{x}^{ + }$ ту часть поверхности $\Gamma $, которая принадлежит плоскости $x = {{\tilde {x}}_{{N + 1}}}$, через $\Gamma _{x}^{ - }$ – часть поверхности $\Gamma $, принадлежащая плоскости $x = {{\tilde {x}}_{0}}$. Аналогично определяются поверхности $\Gamma _{y}^{ + }$, $\Gamma _{y}^{ - }$, $\Gamma _{z}^{ + }$ и $\Gamma _{z}^{ - }$.

Границы отрезка $[a,b]$, на котором идентифицируется коэффициент теплопроводности $K(T)$, задавались как минимальное и максимальное значения функций ${{w}_{0}}(s)$ и ${{w}_{S}}(s,t)$. Отрезок $[a,b]$ разбивался точками ${{\tilde {T}}_{0}} = a$, ${{\tilde {T}}_{1}}$, ${{\tilde {T}}_{2}}$, …, ${{\tilde {T}}_{M}} = b$ на $M$ частей (они могут быть как равными, так и неравными). Каждой из точек ${{\tilde {T}}_{m}}$, $m = 0,\; \ldots ,\;M$, ставилось в соответствие число ${{k}_{m}} = K({{\tilde {T}}_{m}})$. Искомая функция $K(T)$ аппроксимировалась непрерывной кусочно-линейной функцией с узлами в точках $\{ ({{\tilde {T}}_{m}},{{k}_{m}})\} _{{m = 0}}^{M}$ так, как это показано в [17]. Если в процессе решения задачи температура в точке выходила за границы отрезка $[a,b]$, то для определения функции $K(T)$ использовалась линейная экстраполяция.

Итерационная локально-одномерная схема (см. [13]), которая использовалась для решения прямой задачи, является схемой переменных направлений (расщепляется по направлениям $x$, $y$ и $z$). Ниже эта схема приводится в так называемом каноническом виде, который требуется для применения формул быстрого автоматического дифференцирования:

$x$‑направление

$T_{{nil}}^{{j + \frac{1}{3}}} = T_{{nil}}^{j} + {{a}_{{nil}}}\left( {K\left( {T_{{n + 1,il}}^{{j + \frac{1}{3}}}} \right)T_{{n + 1,il}}^{{j + \frac{1}{3}}} + K\left( {T_{{nil}}^{{j + \frac{1}{3}}}} \right)T_{{n + 1,il}}^{{j + \frac{1}{3}}} - K\left( {T_{{n + 1,il}}^{{j + \frac{1}{3}}}} \right)T_{{nil}}^{{j + \frac{1}{3}}} - K\left( {T_{{nil}}^{{j + \frac{1}{3}}}} \right)T_{{nil}}^{{j + \frac{1}{3}}}} \right) + $

$ + \;{{b}_{{nil}}}\left( {K\left( {T_{{n - 1,il}}^{{j + \frac{1}{3}}}} \right)T_{{n - 1,il}}^{{j + \frac{1}{3}}} + K\left( {T_{{nil}}^{{j + \frac{1}{3}}}} \right)T_{{n - 1,il}}^{{j + \frac{1}{3}}} - K\left( {T_{{n - 1,il}}^{{j + \frac{1}{3}}}} \right)T_{{nil}}^{{j + \frac{1}{3}}} - K\left( {T_{{nil}}^{{j + \frac{1}{3}}}} \right)T_{{nil}}^{{j + \frac{1}{3}}}} \right) \equiv \Psi _{{nil}}^{{j + \frac{1}{3}}},$

$y$‑направление

$T_{{nil}}^{{j + \frac{2}{3}}} = T_{{nil}}^{{j + \frac{1}{3}}} + {{c}_{{nil}}}\left( {K\left( {T_{{n,i + 1,l}}^{{j + \frac{2}{3}}}} \right)T_{{n,i + 1,l}}^{{j + \frac{2}{3}}} + K\left( {T_{{nil}}^{{j + \frac{2}{3}}}} \right)T_{{n,i + 1,l}}^{{j + \frac{2}{3}}} - K\left( {T_{{n,i + 1,l}}^{{j + \frac{2}{3}}}} \right)T_{{nil}}^{{j + \frac{2}{3}}} - K\left( {T_{{nil}}^{{j + \frac{2}{3}}}} \right)T_{{nil}}^{{j + \frac{2}{3}}}} \right) + $

$ + \;{{d}_{{nil}}}\left( {K\left( {T_{{n,i - 1,l}}^{{j + \frac{2}{3}}}} \right)T_{{n,i - 1,l}}^{{j + \frac{2}{3}}} + K\left( {T_{{nil}}^{{j + \frac{2}{3}}}} \right)T_{{n,i - 1,l}}^{{j + \frac{2}{3}}} - K\left( {T_{{n,i - 1,l}}^{{j + \frac{2}{3}}}} \right)T_{{nil}}^{{j + \frac{2}{3}}} - K\left( {T_{{nil}}^{{j + \frac{2}{3}}}} \right)T_{{nil}}^{{j + \frac{2}{3}}}} \right) \equiv \Psi _{{nil}}^{{j + \frac{2}{3}}},$

$z$‑направление

$T_{{nil}}^{{j + 1}} = T_{{nil}}^{{j + \frac{2}{3}}} + {{e}_{{nil}}}\left( {K(T_{{ni,l + 1}}^{{j + 1}})T_{{ni,l + 1}}^{{j + 1}} + K(T_{{nil}}^{{j + 1}})T_{{ni,l + 1}}^{{j + 1}} - K(T_{{ni,l + 1}}^{{j + 1}})T_{{nil}}^{{j + 1}} - K(T_{{nil}}^{{j + 1}})T_{{nil}}^{{j + 1}}} \right) + $

$ + \;{{f}_{{nil}}}\left( {K(T_{{ni,l - 1}}^{{j + 1}})T_{{ni,l - 1}}^{{j + 1}} + K(T_{{nil}}^{{j + 1}})T_{{ni,l - 1}}^{{j + 1}} - K(T_{{ni,l - 1}}^{{j + 1}})T_{{nil}}^{{j + 1}} - K(T_{{nil}}^{{j + 1}})T_{{nil}}^{{j + 1}}} \right) \equiv \Psi _{{nil}}^{{j + 1}},$

где

${{a}_{{nil}}} = \frac{\tau }{{6{{C}_{{nil}}}h_{n}^{x}(h_{n}^{x} + h_{{n - 1}}^{x})}},\quad {{b}_{{nil}}} = \frac{\tau }{{6{{C}_{{nil}}}h_{{n - 1}}^{x}(h_{n}^{x} + h_{{n - 1}}^{x})}},\quad {{c}_{{nil}}} = \frac{\tau }{{6{{C}_{{nil}}}h_{i}^{y}(h_{i}^{y} + h_{{i - 1}}^{y})}},$

${{d}_{{nil}}} = \frac{\tau }{{6{{C}_{{nil}}}h_{{i - 1}}^{y}(h_{i}^{y} + h_{{i - 1}}^{y})}},\quad {{e}_{{nil}}} = \frac{\tau }{{6{{C}_{{nil}}}h_{l}^{z}(h_{l}^{z} + h_{{l - 1}}^{z})}},\quad {{f}_{{nil}}} = \frac{\tau }{{6{{C}_{{nil}}}h_{{l - 1}}^{z}(h_{l}^{z} + h_{{l - 1}}^{z})}},$

$(n = \overline {1,\;N - 1,} \quad i = \overline {1,\;I - 1} ,\quad l = \overline {1,\;L - 1} ,\quad j = \overline {0,\;J - 1} ).$

Целевой функционал (1.4) аппроксимировался функцией $F({{k}_{0}},\;{{k}_{1}},\; \ldots ,\;{{k}_{M}})$ конечного числа переменных следующим образом:

(2.1)

$\begin{gathered} \cong \;F = \sum\limits_{j = 1}^J {\left[ {\sum\limits_{l = 1}^{L - 1} {\sum\limits_{i = 1}^{I - 1} {\left( {{{\beta }_{{0il}}}{{{(H_{{0il}}^{j} - P_{{0il}}^{j})}}^{2}} + {{\beta }_{{Nil}}}{{{( - H_{{Nil}}^{j} - P_{{Nil}}^{j})}}^{2}}} \right)h_{i}^{y}h_{l}^{z}} } + } \right.} \\ + \;\sum\limits_{l = 1}^{L - 1} {\sum\limits_{n = 1}^{N - 1} {\left( {{{\beta }_{{n0l}}}{{{(H_{{n0l}}^{j} - P_{{n0l}}^{j})}}^{2}} + {{\beta }_{{nIl}}}{{{( - H_{{nIl}}^{j} - P_{{nIl}}^{j})}}^{2}}} \right)h_{n}^{x}h_{l}^{z}} } + \\ \end{gathered} $

$ + \;\left. {\sum\limits_{n = 1}^{N - 1} {\sum\limits_{i = 1}^{I - 1} {\left( {{{\beta }_{{ni0}}}{{{(H_{{ni0}}^{j} - P_{{ni0}}^{j})}}^{2}} + {{\beta }_{{niL}}}{{{( - H_{{niL}}^{j} - P_{{niL}}^{j})}}^{2}}} \right)h_{n}^{x}h_{i}^{y}} } } \right]\tau .$

Здесь через $H_{{0il}}^{j}$, $ - H_{{Nil}}^{j}$, $H_{{n0l}}^{j}$, $ - H_{{nIl}}^{j}$, $H_{{ni0}}^{j}$ и $ - H_{{niL}}^{j}$ обозначены тепловые потоки, вычисленные на соответствующих гранях параллелепипеда: $\Gamma _{x}^{ - }$, $\Gamma _{x}^{ + }$, $\Gamma _{y}^{ - }$, $\Gamma _{y}^{ + }$, $\Gamma _{z}^{ - }$ и $\Gamma _{z}^{ + }$. Например,

(2.2)

$H_{{0il}}^{j} = \frac{{K(T_{{1il}}^{j}) + K(T_{{0il}}^{j})}}{2}\frac{{T_{{1il}}^{j} - T_{{0il}}^{j}}}{{h_{0}^{x}}},\quad H_{{Nil}}^{j} = \frac{{K(T_{{Nil}}^{j}) + K(T_{{N - 1,il}}^{j})}}{2}\frac{{T_{{Nil}}^{j} - T_{{N - 1,il}}^{j}}}{{h_{{N - 1}}^{x}}}.$

Дискретная сопряженная задача, которая получается в результате применения БАД-методологии в данном случае, совпадает с соотношениями (2.2)–(2.5), которые представлены в [10], однако формулы для вычисления производных $\frac{{\partial F}}{{\partial T_{{nil}}^{j}}}$ будут теперь иметь такой вид:

$\frac{{\partial F}}{{\partial T_{{nil}}^{j}}} = 0,\quad n = \overline {2,\;N - 2} ,\quad i = \overline {2,\;I - 2} ,\quad l = \overline {2,\;L - 2} ,$

$\frac{{\partial F}}{{\partial T_{{0il}}^{j}}} = \tilde {A}_{{il}}^{j}\left[ {K{\text{'}}(T_{{0il}}^{j})(T_{{1il}}^{j} - T_{{0il}}^{j}) - (K(T_{{0il}}^{j}) + K(T_{{1il}}^{j}))} \right],\quad i = \overline {1,\;I - 1} ,\quad l = \overline {1,\;L - 1} ,$

$\frac{{\partial F}}{{\partial T_{{1il}}^{j}}} = \tilde {A}_{{il}}^{j}\left[ {K{\text{'}}(T_{{1il}}^{j})(T_{{1il}}^{j} - T_{{0il}}^{j}) + (K(T_{{0il}}^{j}) + K(T_{{1il}}^{j}))} \right],\quad i = \overline {2,\;I - 1} ,\quad l = \overline {2,\;L - 1} ,$

$\frac{{\partial F}}{{\partial T_{{N - 1,il}}^{j}}} = \tilde {B}_{{il}}^{j}\left[ {K{\text{'}}(T_{{N - 1,il}}^{j})(T_{{Nil}}^{j} - T_{{N - 1,il}}^{j}) - (K(T_{{Nil}}^{j}) + K(T_{{N - 1,il}}^{j}))} \right],\quad i = \overline {1,\;I - 2} ,\quad l = \overline {1,\;L - 2} ,$

$\frac{{\partial F}}{{\partial T_{{Nil}}^{j}}} = \tilde {B}_{{il}}^{j}\left[ {K{\text{'}}(T_{{Nil}}^{j})(T_{{Nil}}^{j} - T_{{N - 1,il}}^{j}) + (K(T_{{Nil}}^{j}) + K(T_{{N - 1,il}}^{j}))} \right],\quad i = \overline {1,\;I - 1} ,\quad l = \overline {1,\;L - 1} ,$

$\tilde {A}_{{il}}^{j} = \frac{{{{\beta }_{{0il}}}h_{i}^{y}h_{l}^{z}\tau }}{{h_{0}^{x}}}(H_{{0il}}^{j} - P_{{0il}}^{j}),\quad \tilde {B}_{{il}}^{j} = \frac{{{{\beta }_{{Nil}}}h_{i}^{y}h_{l}^{z}\tau }}{{h_{{N - 1}}^{x}}}(H_{{Nil}}^{j} + P_{{Nil}}^{j}),$

$K{\text{'}}(T_{{nil}}^{j}) = \frac{{dK(T)}}{{dT}}(T_{{nil}}^{j}) = \frac{{{{k}_{m}} - {{k}_{{m - 1}}}}}{{{{{\tilde {T}}}_{m}} - {{{\tilde {T}}}_{{m - 1}}}}},$

а индекс $m$ определяется условием ${{\tilde {T}}_{{m - 1}}} \leqslant T \leqslant {{\tilde {T}}_{m}}$.

Производные

$\frac{{\partial F}}{{\partial T_{{n0l}}^{j}}},\quad n = \overline {1,N - 1} ,\quad l = \overline {1,L - 1} ,\quad \frac{{\partial F}}{{\partial T_{{n1l}}^{j}}},\quad n = \overline {2,N - 1} ,\quad l = \overline {2,L - 1} ,$

$\frac{{\partial F}}{{\partial T_{{n,I - 1,l}}^{j}}},\quad n = \overline {1,\;N - 2} ,\quad l = \overline {1,\;L - 2} ,\quad \frac{{\partial F}}{{\partial T_{{nIl}}^{j}}},\quad n = \overline {1,\;N - 1} ,\quad l = \overline {1,\;L - 1} ,$

$\frac{{\partial F}}{{\partial T_{{ni0}}^{j}}},\quad n = \overline {1,\;N - 1} ,\quad i = \overline {1,\;I - 1} ,\quad \frac{{\partial F}}{{\partial T_{{ni1}}^{j}}},\quad n = \overline {2,\;N - 1} ,\quad i = \overline {2,\;I - 1} ,$

$\frac{{\partial F}}{{\partial T_{{ni,L - 1}}^{j}}},\quad n = \overline {1,N - 2} ,\quad i = \overline {1,I - 2} \quad {\text{и}}\quad \frac{{\partial F}}{{\partial T_{{niL}}^{j}}},\quad n = \overline {1,\;N - 1} ,\quad i = \overline {1,\;I - 1} $

вычисляются подобным образом.

При вычислении производных, в которых встречаются повторяющиеся пространственные индексы, такие как:

$\frac{{\partial F}}{{\partial T_{{11l}}^{j}}},\quad l = \overline {2,\;L - 1} ,\quad \frac{{\partial F}}{{\partial T_{{n11}}^{j}}},\quad n = \overline {2,\;N - 1} ,$

$\frac{{\partial F}}{{\partial T_{{1i1}}^{j}}},\quad i = \overline {2,\;I - 1} ,\quad \frac{{\partial F}}{{\partial T_{{N - 1,I - 1,l}}^{j}}},\quad l = \overline {1,\;L - 1} \quad {\text{и}}\quad N = I,$

$\frac{{\partial F}}{{\partial T_{{n,I - 1,L - 1}}^{j}}},\quad n = \overline {1,\;N - 1} \quad {\text{и}}\quad I = L,\quad \frac{{\partial F}}{{\partial T_{{N - 1,i,L - 1}}^{j}}},\quad i = \overline {1,\;I - 1} \quad {\text{и}}\quad N = L,$

необходимо учитывать все пересечения по этим индексам. Например,

$\begin{gathered} \frac{{\partial F}}{{\partial T_{{11l}}^{j}}} = \tilde {A}_{{1l}}^{j}\left[ {K{\text{'}}(T_{{11l}}^{j})(T_{{11l}}^{j} - T_{{01l}}^{j}) + (K(T_{{01l}}^{j}) + K(T_{{11l}}^{j}))} \right] + \\ + \;\tilde {C}_{{1l}}^{j}\left[ {K{\text{'}}(T_{{11l}}^{j})(T_{{11l}}^{j} - T_{{10l}}^{j}) + (K(T_{{10l}}^{j}) + K(T_{{11l}}^{j}))} \right],\quad l = \overline {2,\;L - 1} , \\ \end{gathered} $

$\begin{gathered} \frac{{\partial F}}{{\partial T_{{111}}^{j}}} = \tilde {A}_{{11}}^{j}\left[ {K{\text{'}}(T_{{111}}^{j})(T_{{111}}^{j} - T_{{011}}^{j}) + (K(T_{{011}}^{j}) + K(T_{{111}}^{j}))} \right] + \\ + \;\tilde {C}_{{11}}^{j}\left[ {K{\text{'}}(T_{{111}}^{j})(T_{{111}}^{j} - T_{{101}}^{j}) + (K(T_{{101}}^{j}) + K(T_{{111}}^{j}))} \right] + \\ \end{gathered} $

$\begin{gathered} + \;\tilde {E}_{{11}}^{j}\left[ {K{\text{'}}(T_{{111}}^{j})(T_{{111}}^{j} - T_{{110}}^{j}) + (K(T_{{110}}^{j}) + K(T_{{111}}^{j}))} \right], \\ \tilde {C}_{{nl}}^{j} = \frac{{{{\beta }_{{n0l}}}h_{n}^{x}h_{l}^{z}\tau }}{{h_{0}^{y}}}(H_{{n0l}}^{j} - P_{{n0l}}^{j}),\quad \tilde {E}_{{ni}}^{j} = \frac{{{{\beta }_{{ni0}}}h_{n}^{x}h_{i}^{y}\tau }}{{h_{0}^{z}}}(H_{{ni0}}^{j} - P_{{ni0}}^{j}). \\ \end{gathered} $

Компоненты градиента целевой функции, в соответствии с БАД-методологией, вычисляются по формуле:

$\frac{{\partial F}}{{\partial {{k}_{m}}}} = \sum\limits_{g = 0}^J {\sum\limits_{r = 0}^L {\sum\limits_{q = 0}^I {\sum\limits_{s = 0}^N {\left( {{\rm X}_{{rqs}}^{g}\frac{{\partial K(T_{{rqs}}^{g})}}{{\partial {{k}_{m}}}}} \right)} } } } + \sum\limits_{g = 0}^{J - 1} {\sum\limits_{r = 0}^L {\sum\limits_{q = 0}^I {\sum\limits_{s = 0}^N {} } } } \left( {{\rm X}_{{rqs}}^{{g + \frac{1}{3}}}\frac{{\partial K\left( {T_{{rqs}}^{{g + \frac{1}{3}}}} \right)}}{{\partial {{k}_{m}}}} + {\rm X}_{{rqs}}^{{g + \frac{2}{3}}}\frac{{\partial K\left( {T_{{rqs}}^{{g + \frac{2}{3}}}} \right)}}{{\partial {{k}_{m}}}}} \right) + $

(2.3)

$\begin{gathered} + \;\sum\limits_{j = 1}^J {\left[ {\sum\limits_{l = 1}^{L - 1} {\sum\limits_{i = 1}^{I - 1} {\left( {\frac{{\partial F}}{{\partial K(T_{{0il}}^{j})}} + \frac{{\partial F}}{{\partial K(T_{{1il}}^{j})}} + \frac{{\partial F}}{{\partial K(T_{{N - 1,il}}^{j})}} + \frac{{\partial F}}{{\partial K(T_{{Nil}}^{j})}}} \right)} } + } \right.} \\ + \;\sum\limits_{l = 1}^{L - 1} {\sum\limits_{n = 1}^{N - 1} {\left( {\frac{{\partial F}}{{\partial K(T_{{n0l}}^{j})}} + \frac{{\partial F}}{{\partial K(T_{{n1l}}^{j})}} + \frac{{\partial F}}{{\partial K(T_{{n,I - 1,l}}^{j})}} + \frac{{\partial F}}{{\partial K(T_{{nIl}}^{j})}}} \right)} } + \\ \end{gathered} $

$ + \;\left. {\sum\limits_{n = 1}^{N - 1} {\sum\limits_{i = 1}^{I - 1} {\left( {\frac{{\partial F}}{{\partial K(T_{{ni0}}^{j})}} + \frac{{\partial F}}{{\partial K(T_{{ni1}}^{j})}} + \frac{{\partial F}}{{\partial K(T_{{ni,L - 1}}^{j})}} + \frac{{\partial F}}{{\partial K(T_{{niL}}^{j})}}} \right)} } } \right],\quad m = \overline {0,\;M} .$

Производные $\frac{{\partial F}}{{\partial K(T_{{nil}}^{j})}}$ вычисляются по формулам:

$\frac{{\partial F}}{{\partial K(T_{{0il}}^{j})}} = \frac{{\partial F}}{{\partial K(T_{{1il}}^{j})}} = \tilde {A}_{{il}}^{j}(T_{{1il}}^{j} - T_{{0il}}^{j}),\quad \frac{{\partial F}}{{\partial K(T_{{Nil}}^{j})}} = \frac{{\partial F}}{{\partial K(T_{{N - 1,il}}^{j})}} = \tilde {B}_{{il}}^{j}(T_{{Nil}}^{j} - T_{{N - 1,il}}^{j}).$

Остальные подобные производные вычисляются аналогичным образом.

Выражения для ${\rm X}_{{rqs}}^{g}$, ${\rm X}_{{rqs}}^{{g + 1/3}}$, ${\rm X}_{{rqs}}^{{g + 2/3}}$, содержащие сопряженные переменные, приведены в работе [10], а множители $\frac{{\partial K(T_{{rqs}}^{g})}}{{\partial {{k}_{m}}}}$, присутствующие в формуле (2.3), должны вычисляться так:

$\frac{{\partial K(T_{{ni}}^{j})}}{{\partial {{k}_{{m - 1}}}}} = 1 - \frac{{T_{{ni}}^{j} - {{{\tilde {T}}}_{{m - 1}}}}}{{{{{\tilde {T}}}_{m}} - {{{\tilde {T}}}_{{m - 1}}}}},\quad \frac{{\partial K(T_{{ni}}^{j})}}{{\partial {{k}_{m}}}} = \frac{{T_{{ni}}^{j} - {{{\tilde {T}}}_{{m - 1}}}}}{{{{{\tilde {T}}}_{m}} - {{{\tilde {T}}}_{{m - 1}}}}},$

где индекс $m$ определяется условием ${{\tilde {T}}_{{m - 1}}} \leqslant T \leqslant {{\tilde {T}}_{m}}$, $m = 1,\; \ldots ,\;M$.

Отметим, что значение градиента целевой функции $F({{k}_{0}},\;{{k}_{1}},\; \ldots ,\;{{k}_{M}})$, вычисленное по формуле (2.3), является точным для выбранной аппроксимации задачи оптимального управления, что является чрезвычайно важным при использовании градиентных методов.

3. РЕЗУЛЬТАТЫ ЧИСЛЕННЫХ РАСЧЕТОВ

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

Для получения “экспериментального” потока тепла $P(x,y,z,t)$ на поверхности параллелепипеда решалась прямая начально-краевая задача (1.1)–(1.3) с заданным коэффициентом теплопроводности, определялось поле температур в объекте, а затем по нему рассчитывались потоки с помощью формул, аналогичных формулам (2.2).

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

Выбранная в настоящей работе локально-одномерная схема аппроксимации уравнения теплопроводности является устойчивой и позволяет работать с достаточно большим шагом по времени. Тем не менее для каждой используемой пространственной сетки необходимо проводить исследования, касающиеся выбора сетки по времени. Представленные здесь результаты численных расчетов получены при использовании равномерной пространственной сетки с параметрами $N = I = L = 25$. Количество интервалов разбиения отрезка $[0,\Theta ]$ по времени $J = 25$.

Минимизация целевой функции $F({{k}_{0}},\;{{k}_{1}},\; \ldots ,\;{{k}_{M}})$ проводилась численно с помощью градиентного метода. Отрезок $[a,b]$, на котором идентифицировался коэффициент теплопроводности, разбивался на 64 интервала, т.е. $M = 64$. Хочется отметить, что для ускорения итерационного процесса спуска по градиенту, как и в работе [10], во всех представленных здесь примерах использовался подход, предложенный в [6]. Он основан на последовательном увеличении числа $M$ разбиений отрезка $[a,b]$. Начинать процесс желательно с $M = 1$. После получения оптимального решения его следует использовать в качестве начального приближения для варианта с $M = 2$. Получив оптимальное решение для $M = 2$, использовать его в качестве начального приближения для $M = 4$, и т.д.

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

(3.1)

${{\varepsilon }_{1}} = \mathop {\max }\limits_{0 \leqslant m \leqslant M} \frac{{\left| {{{K}_{{{\text{opt}}}}}({{{\tilde {T}}}_{m}}) - K({{{\tilde {T}}}_{m}})} \right|}}{{K{\text{*}}}}$

(3.2)

${{\varepsilon }_{2}} = \frac{1}{{K{\text{*}}}}\sqrt {\sum\limits_{m = 0}^M {\frac{{{{{({{K}_{{{\text{opt}}}}}({{{\tilde {T}}}_{m}}) - K({{{\tilde {T}}}_{m}}))}}^{2}}}}{{M + 1}}} } ,$

где $K({{\tilde {T}}_{m}})$ – значения аналитического коэффициента теплопроводности, вычисленные в опорных точках отрезка температур, ${{K}_{{{\text{opt}}}}}({{\tilde {T}}_{m}})$ – значения коэффициента теплопроводности, полученные в результате решения оптимизационной задачи, а $K{\text{*}}$ – некоторое характерное значение функции $K(T)$. В работе полагалось, что

$K{\kern 1pt} * = \frac{{\sum\nolimits_{m = 0}^M {K({{{\tilde {T}}}_{m}})} }}{{M + 1}}$.

Проведенные численные эксперименты показали, что качество восстановления коэффициента теплопроводности сильно зависит от плотности распределения “экспериментальных” данных. Бывают случаи, когда на некоторых участках отрезка $[a,b]$ слишком мало данных, необходимых для идентификации коэффициента теплопроводности. Поэтому в каждом конкретном случае необходимо проанализировать распределение “экспериментальных” данных по интервалам интересующего нас отрезка температур.

Для анализа распределения экспериментальных данных по интервалам интересующего нас отрезка температур введем функцию ${{W}_{*}}(T)$ – относительную меру той подобласти области $Q \times (0,\Theta )$, температура в точках которой меньше $T$, а через $W(T)$ обозначим плотность (производную) функции ${{W}_{*}}(T)$ по переменной $T$.

Для надежности все рассмотренные задачи решались численно несколько раз, при этом в качестве начального приближения ${{K}_{{ini}}}(T)$ выбирались разные функции.

3.1. Первая серия расчетов

В первой серии расчетов рассматривалась задача нахождения коэффициента теплопроводности вещества при следующих входных данных: в качестве начальной функции ${{w}_{0}}(x,y,z)$ и граничной функции ${{w}_{\Gamma }}(x,y,z,t)$ выбирались следы функции

(3.3)

$\Lambda (x,y,z,t) = x + y + z + 3t + 0.5$

на параболической границе области $Q \times (0,\Theta ) = (0,\;1) \times (0,\;1) \times (0,\;1) \times (0,\;1)$. В этом случае функция (3.3) является решением смешанной задачи (1.1)–(1.3) при $C(s) = 1$ и $K(T) = T$. Температура на параболической границе рассматриваемой области изменяется от $a = {\text{0}}{\text{.5}}$ до $b = 6.{\text{5}}$.

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

Все проведенные расчеты первой серии можно разделить на несколько групп. В расчетах первой группы полагалось, что в целевом функционале весовой параметр $\beta (s(\Gamma ))$ равняется $1$ во всех точках границы $\Gamma $ объекта. На фиг. 1 представлено распределение функции $W(T)$ вдоль отрезка $[0.5,\;6.5]$. Оно напоминает нормальное распределение с центром в точке $T = 3.5$. Отрезку температур $T \in [2.5,\;5.0]$ соответствует основное количество экспериментальных данных, а в окрестности концов отрезка $[0.5,\;6.5]$ этих данных сравнительно мало.

Фиг. 1.

Проведенные расчеты показали, что предложенный алгоритм приводит к одному и тому же оптимальному решению независимо от того, какая функция выбрана в качестве начального приближения. Во всех случаях коэффициент теплопроводности идентифицировался с машинной точностью. В процессе решения оптимизационной задачи целевой функционал изменился от ${\text{2}}{\text{.4159}} \times {{10}^{1}}$ на начальном приближении (${{K}_{{ini}}}(T) = 5.0$) к коэффициенту теплопроводности до ${\text{1}}{\text{.5607}} \times {{10}^{{ - 24}}}$ на финальном (оптимальном) приближении, при этом максимум модуля градиента уменьшился от ${\text{1}}{\text{.6057}}$ до ${\text{9}}{\text{.4125}} \times {{10}^{{ - 14}}}$.

Отметим, что когда в качестве “экспериментальных” данных использовалось температурное поле (см. [10]), то обратная задача, рассмотренная в этой серии, имела неединственное решение и для получения единственного решения требовалось задание дополнительного условия. Отметим также, что неравномерное распределение экспериментальных данных сказалось здесь только на сходимости итерационного процесса, но не на качестве полученного решения.

К расчетам второй группы относятся те случаи, когда в качестве “экспериментального” теплового потока использовался поток тепла, вычисленный только на одной из шести граней параллелепипеда. На фиг. 2 представлено распределение функции $W(T)$ вдоль отрезка $[0.5,\;6.5]$ для случая, когда в качестве “экспериментальных” данных использовался поток тепла только на той грани, для которой $z = 0$, т.е. ${{\beta }_{{ni0}}} = 1$ для всех $n = \overline {1,\;N - 1} $, $i = \overline {1,\;I - 1} $ (на остальных гранях весовая функция равна нулю). Качественно функция $W(T)$ здесь похожа на ту, что представлена на фиг. 1, но на правом конце отрезка вообще отсутствуют экспериментальные данные. Тем не менее предложенный алгоритм решения обратной задачи позволил восстановить коэффициент теплопроводности с машинной точностью на всем отрезке $[0.5,\;6.5]$. По-видимому, это связано с линейным характером искомой функции $K(T) = T$.

Фиг. 2.

Распределение функции $W(T)$, подобное приведенному на фиг. 2, получается и тогда, когда в качестве “экспериментальных” данных используются потоки тепла на тех гранях параллелепипеда, для которых $x = 0$ или $y = 0$. И в случаях использования экспериментальных данных только на этих гранях по отдельности удается восстановить коэффициент теплопроводности с машинной точностью.

3.2. Вторая серия раcчетов

Во второй серии расчетов рассматривалась задача определения коэффициента теплопроводности при таких входных данных: в качестве начальной функции ${{w}_{0}}(x,y,z)$ и граничной функции ${{w}_{\Gamma }}(x,y,z,t)$ выбирались следы функции

(3.4)

$\Lambda (x,y,z,t) = \sqrt {\frac{{{{x}^{2}} + {{y}^{2}} + {{z}^{2}}}}{{9 - 8t}}} $

на параболической границе области $Q \times (0,\Theta ) = (0,\;1) \times (0,\;1) \times (0,\;1) \times (0,\;1)$. Функция (3.4) является в этом случае решением смешанной задачи (1.1)–(1.3) при $C(s) = 1$ и $K(T) = {{T}^{2}}$. Температура на параболической границе рассматриваемой области изменяется от $a = 0.0$ до $b = 1.732$.

В расчетах первой группы полагалось, что в целевом функционале весовой параметр $\beta (s(\Gamma ))$ равняется $1$ во всех точках границы $\Gamma $ объекта. Распределение экспериментальных данных вдоль отрезка температур $[0.0,\;1.732]$ в этом случае представлено на фиг. 3 в виде зависимости функции $W(T)$ от температуры. Анализ этого распределения показывает, что правому концу отрезка $[0.0,\;1.732]$ соответствует слишком мало “экспериментальных” данных. Поэтому можно предположить, что там возникнут трудности с восстановлением коэффициента теплопроводности. Проведенные расчеты подтвердили это предположение: при равномерном разбиении отрезка температур на $M = 64$ интервалов последняя компонента градиента функционала не изменяется (она равна нулю с машинной точностью).

Фиг. 3.

Процесс построения решения задачи идентификации в том случае, когда в качестве начального приближения было выбрано управление ${{K}_{{ini}}}(T) \equiv 5.0$, проиллюстрирован на фиг. 4. Здесь представлены функция $K(T) = {{T}^{2}}$ (теоретическое значение коэффициента теплопроводности) и оптимальные управления, полученные при $M = 1,\;2,\;4$. Видно, что опорные точки кусочно-линейного оптимального управления, полученного при $M = 4$, почти лежат на линии $K(T) = {{T}^{2}}$. Что касается оптимального управления, полученного при $M = 8$, то оно практически совпадает с функцией $K(T) = {{T}^{2}}$ (${{\varepsilon }_{1}} = {\text{9}}{\text{.4417}} \times {\text{1}}{{{\text{0}}}^{{ - {\text{3}}}}}$ и ${{\varepsilon }_{2}} = {\text{7}}{\text{.1752}} \times {\text{1}}{{{\text{0}}}^{{ - {\text{3}}}}}$). Отклонения полученного при $M = 64$ коэффициента теплопроводности ${{K}_{{{\text{opt}}}}}(T)$ от его аналитического значения $K(T) = {{T}^{2}}$ составляют ${{\varepsilon }_{1}} = {\text{5}}{\text{.3197}} \times {\text{1}}{{{\text{0}}}^{{ - {\text{4}}}}}$ и ${{\varepsilon }_{2}} = {\text{7}}{\text{.2643}} \times {\text{1}}{{{\text{0}}}^{{ - {\text{5}}}}}$ соответственно.

Фиг. 4.

К расчетам второй группы относятся те случаи, когда в качестве “экспериментального” теплового потока использовался поток тепла, вычисленный на двух гранях параллелепипеда: при $x = 0$ и $x = 1$, т.е. в (2.1) ${{\beta }_{{0il}}} = 1$ и ${{\beta }_{{Nil}}} = 1$ для всех $i = \overline {1,\;I - 1} $, $l = \overline {1,\;L - 1} $, а в остальных точках весовая функция равна нулю.

Распределение функции $W(T)$ вдоль отрезка $[0.0,\;1.732]$ близко к тому, которое представлено на фиг. 3. Поэтому можно ожидать, что в этом случае коэффициент теплопроводности будет восстанавливаться с той же точностью, которая была получена при использовании “экспериментального” теплового потока на всех шести гранях параллелепипеда. Проведенные численные расчеты подтвердили это предположение. Нельзя не отметить, что в этом случае оптимизационная задача решается значительно быстрее.

3.3. Третья серия расчетов

В третьей серии расчетов задача нахождения коэффициента теплопроводности решалась при следующих входных данных: в качестве начальной функции ${{w}_{0}}(x,y,z)$ и граничной функции ${{w}_{\Gamma }}(x,y,z,t)$ выбирались следы функции

(3.5)

$\Lambda (x,y,z,t) = \frac{3}{{1.8(5 - x - y - z - 1.8t)}}$

на параболической границе области $Q \times (0,\Theta ) = (0,\;1) \times (0,\;1) \times (0,\;1) \times (0,\;1)$. Функция (3.5) является решением смешанной задачи (1.1)–(1.3) при $C(s) = 1$ и $K(T) = \frac{1}{T}$. Температура на параболической границе рассматриваемой области изменяется от $a = {\text{0}}{\text{.33}}$ до $b = 8.{\text{33}}$.

В расчетах первой группы этой серии полагалось, что в целевом функционале весовой параметр $\beta (s(\Gamma ))$ равняется $1$ во всех точках границы $\Gamma $ объекта. На фиг. 5 представлено распределение “экспериментальных” данных по отрезку температур $[0.33,\;8.33]$ (функция $W(T)$). Анализируя представленный на этой фигуре график, нетрудно увидеть, что на правом конце отрезка $[0.33,\;8.33]$ практически нет “экспериментальных” данных. Как показали проведенные расчеты этой серии, коэффициент теплопроводности на подотрезке $[6.0,\;8.33]$ не восстанавливается. При разбиении отрезка температур на $M = 64$ интервала последние 19 компонент градиента функционала всегда оказывались нулевыми (с машинной точностью).

Фиг. 5.

Задача оптимизации, к которой сводится рассматриваемая обратная коэффициентная задача, решалась при разных начальных приближениях ${{K}_{{ini}}}(T)$. Процесс построения решения задачи идентификации в том случае, когда в качестве начального приближения было выбрано управление ${{K}_{{ini}}}(T) \equiv 5.0$, проиллюстрирован на фиг. 6. Здесь представлены функция $K(T) = \frac{1}{T}$ (теоретическое значение коэффициента теплопроводности) и оптимальные управления, полученные при $M = 2,\;4,\;16$. График оптимального управления, полученного при $M = 64$ на отрезке $[0.33,\;6.0]$, полностью совпадает с графиком функции $K(T) = \frac{1}{T}$. Отклонения полученного коэффициента теплопроводности ${{K}_{{{\text{opt}}}}}(T)$ от его аналитического значения $K(T) = \frac{1}{T}$ на отрезке $[0.33,\;6.0]$ составляют ${{\varepsilon }_{1}} = 3.{\text{5506}} \times {\text{1}}{{{\text{0}}}^{{ - {\text{4}}}}}$ и ${{\varepsilon }_{2}} = 1.{\text{5300}} \times {\text{1}}{{{\text{0}}}^{{ - {\text{4}}}}}$ соответственно. Однако на отрезке $[6.0,\;8.33]$ сохранился след от управления, полученного при $M = 16$, не совпадающий с $K(T) = \frac{1}{T}$.

Фиг. 6.

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

К расчетам второй группы относятся те случаи, когда в качестве “экспериментального” теплового потока использовался поток тепла, вычисленный только на одной грани параллелепипеда, а именно на той, для которой $x = 0$, т.е. в (2.1) ${{\beta }_{{0il}}} = 1$ для всех $i = \overline {1,\;I - 1} $, $l = \overline {1,\;L - 1} $, а в остальных точках весовая функция равна нулю. На фиг. 7 представлено распределение функции $W(T)$ вдоль отрезка $[0.33,\;8.33]$. В этом случае на отрезке температур $[1.5,\;8.33]$ практически нет “экспериментальных” данных.

Фиг. 7.

Построение решения задачи идентификации проводилось обычным образом. В качестве начального приближения было выбрано управление ${{K}_{{ini}}}(T) \equiv 5.0$. Затем задача решалась последовательно при $M = 1,\;2,\; \ldots ,\;64$. График оптимального управления, полученного при $M = 64$, полностью совпадает с графиком функции $K(T) = \frac{1}{T}$ только на отрезке $[0.33,\;1.5]$ (см. фиг. 8), т.е. для тех значений температур, при которых задан тепловой поток. Отклонения полученного коэффициента теплопроводности ${{K}_{{{\text{opt}}}}}(T)$ от его аналитического значения $K(T) = \frac{1}{T}$ на отрезке $[0.33,\;1.5]$ составляют ${{\varepsilon }_{1}} = {\text{5}}{\text{.8421}} \times {\text{1}}{{{\text{0}}}^{{ - {\text{4}}}}}$ и ${{\varepsilon }_{2}} = {\text{2}}{\text{.0151}} \times {\text{1}}{{{\text{0}}}^{{ - {\text{4}}}}}$ соответственно. На отрезке $[0.33,\;8.33]$ сохранился след управления, полученного при $M = 16$.

Фиг. 8.

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

3.4. Четвертая серия расчетов

В четвертой серии расчетов рассматривалась задача нахождения коэффициента теплопроводности вещества при следующих входных данных: в качестве начальной функции ${{w}_{0}}(x,y,z)$ и граничной функции ${{w}_{\Gamma }}(x,y,z,t)$ выбирались следы функции

$\Lambda (x,y,z,t) = \sqrt {\frac{{9{{{(x + 1)}}^{2}} + 20{{y}^{2}} + 25{{z}^{2}}}}{{9 - 8t}}} $

на параболической границе области $Q \times (0,\Theta ) = (0,\;1) \times (0,\;1) \times (0,\;1) \times (0,\;1)$.

В качестве “экспериментальных” данных использовался поток тепла, вычисленный по температурному полю, полученному в результате решения прямой задачи (1.1)–(1.3) при $C(s) = 1$ и с коэффициентом теплопроводности $K(T) = k(T)$, где функция $k(T)$ определялась следующим равенством:

$k(T) = \left\{ \begin{gathered} 0.1(T - 3)(T - 6)(T - 7) + 3.4,\quad T \geqslant 3, \hfill \\ 1.2(T - 3) + 3.4,\quad T < 3. \hfill \\ \end{gathered} \right.$

Анализ построенного таким образом температурного поля позволил определить диапазон изменения температуры: $a = 1$, $b = 9$.

В первом примере этой серии расчетов полагалось, что в целевом функционале весовой параметр $\beta (s(\Gamma ))$ равняется $1$ во всех точках границы $\Gamma $ объекта. На фиг. 9 представлено распределение функции $W(T)$ вдоль отрезка $[1.0,\;9.0]$. В данном случае практически нет “экспериментальных” данных лишь на небольшом подотрезке на правом конце отрезка температур $[1.0,\;9.0]$. Можно предположить, что здесь возникнут трудности с восстановлением коэффициента теплопроводности. Проведенные расчеты подтвердили это предположение: при равномерном разбиении отрезка температур на $M = 64$ интервала последняя компонента градиента функционала не изменялась (она равнялась нулю с машинной точностью).

Фиг. 9.

Процесс построения решения задачи идентификации в том случае, когда в качестве начального приближения было выбрано управление ${{K}_{{ini}}}(T) \equiv 5.0$, проиллюстрирован на фиг. 10. Здесь представлены функция $k(T)$ и оптимальные управления, полученные при $M = 2,\;4,\;8$. Видно, что опорные точки кусочно-линейного оптимального управления, полученного при $M = 8$, почти лежат на линии $K(T) = k(T)$.

Фиг. 10.

График оптимального управления, полученного при $M = 64$, совпадает с графиком функции $k(T)$. Отклонения полученного при $M = 64$ коэффициента теплопроводности ${{K}_{{{\text{opt}}}}}(T)$ от его аналитического значения $k(T)$ на отрезке $[1.0,\;8.875]$ составляют ${{\varepsilon }_{1}} = {\text{1}}{\text{.3914}} \times {\text{1}}{{{\text{0}}}^{{ - {\text{7}}}}}$ и ${{\varepsilon }_{2}} = {\text{1}}{\text{.8490}} \times {\text{1}}{{{\text{0}}}^{{ - {\text{8}}}}}$ соответственно. На последнем интервале сохраняется след значения коэффициента теплопроводности, полученного при последнем разбиении отрезка температур, для которого соответствующая компонента градиента целевого функционала ненулевая.

Во втором примере четвертой серии расчетов полагалось, что в качестве “экспериментального” теплового потока использовался поток тепла, вычисленный на двух гранях параллелепипеда: при $x = 0$ и $x = 1$, т.е. в (2.1) ${{\beta }_{{0il}}} = 1$ и ${{\beta }_{{Nil}}} = 1$ для всех $i = \overline {1,\;I - 1} $, $l = \overline {1,\;L - 1} $, а в остальных точках весовая функция равна нулю.

Распределение функции $W(T)$ вдоль отрезка $[1.0,\;9.0]$ практически совпадает с тем, которое представлено на фиг. 9. В этом случае можно ожидать, что коэффициент теплопроводности будет восстанавливаться с той же точностью, которая была получена при использовании “экспериментального” теплового потока на всех шести гранях параллелепипеда. Проведенные численные расчеты подтвердили это предположение, а также тот факт, что оптимизационная задача решается быстрее, если использовать тепловые потоки только на этих двух гранях.

Полученные результаты позволяют использовать улучшенную стратегию при решении обратной коэффициентной задачи. Вначале надо построить распределение функции $W(T)$ на основе экспериментальных данных на всех гранях параллелепипеда, а затем на каждой грани по отдельности. Если окажется, что распределение функции $W(T)$ на какой-то грани качественно похоже на ее распределение по всей поверхности, то при решении обратной задачи достаточно использовать экспериментальные данные только на этой одной грани. Поиск решения обратной задачи при этом значительно упростится.

Как и при решении обратной коэффициентной задачи меньшей размерности, предложенный алгоритм продемонстрировал свою устойчивость: небольшие отклонения в экспериментальных данных (~10%) приводят к погрешностям в решении того же порядка (∼2–5%). При бóльших отклонениях в экспериментальных данных для получения устойчивого решения необходимо вводить регуляризатор (см. [9]).

ЗАКЛЮЧЕНИЕ

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

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

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

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

Особо хочется остановиться на сравнении результатов решения задачи идентификации коэффициента теплопроводности в двух различных постановках: использование в качестве “экспериментальных” данных заданного температурного поля (функционал “поле”) или заданного теплового потока на границе объекта (функционал “поток”). Как следует из работы [10], использование функционала “поле” может приводить к неединственности решения задачи. При исследовании большого числа задач с функционалом “поток” никогда не приходилось сталкиваться с неединственностью решения. В пользу использования функционала “поток” при решении задач идентификации коэффициента теплопроводности говорит и то, что тепловой поток на границе объекта замерить проще, чем температуру в самом объекте.

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

Алифанов О.М., Черепанов В.В. Математическое моделирование высокопористых волокнистых материалов и определение их физических свойств // Теплофиз. высоких температур. 2009. Т. 47. № 3. С. 463–472.
Алифанов О.М. Обратные задачи теплообмена. М.: Машиностр., 1988.
Зубов В.И. Применение методологии быстрого автоматического дифференцирования к решению обратной коэффициентной задачи для уравнения теплопроводности // Ж. вычисл. матем. и матем. физ. 2016. Т. 56. № 10. С. 1760–1774.
Albu A.F., Evtushenko Y.G., Zubov V.I. Identification of Discontinuous Thermal Conductivity Coefficient Using Fast Automatic Differentiation. In: Battiti R., Kvasov D., Sergeyev Y. (eds) Learning and Intelligent Optimization. LION 2017. Lecture Notes in Computer Science, 2017. V. 10556. P. 295–300. Springer, Cham.
Zubov V.I., Albu A.F. The FAD-methodology and Recovery the Thermal Conductivity Coefficient in Two Dimension Case. Proc. of the VIII Internat. Conference on Optimizat. Methods and Applications “Optimization and applications”, 2017. P. 39–44.
Албу А.Ф., Зубов В.И. О восстановлении коэффициента теплопроводности вещества по температурному полю // Ж. вычисл. матем. и матем. физ. 2018. Т. 58. № 10. С. 1642–1657.
Албу А.Ф., Зубов В.И. Восстановление коэффициента теплопроводности вещества по поверхностному тепловому потоку // Ж. вычисл. матем. и матем. физ. 2018. Т. 58. № 12. С. 1992–2017.
Alla Albu, Vladimir Zubov. Identification of the thermal conductivity coefficient in two dimension case // Optim. Lett. 2018. P. 1–17.
Albu A., Zubov V. On the Stability of the Algorithm of Identification of the Thermal Conductivity Coefficient. In: Evtushenko Y., Jaćimović M., Khachay M., Kochetov Y., Malkova V., Posypkin M. (Eds.) Optimization and Applications. OPTIMA 2018. Communications in Computer and Information Science. 2019. V. 974. Springer, Cham. P. 247–263.
Албу А.Ф., Зубов В.И. Идентификация коэффициента теплопроводности вещества в трехмерном случае путем решения соответствующей задачи оптимизации // Ж. вычисл. матем. и матем. физ. 2021. Т. 61. № 9. С. 1447–1463.
Евтушенко Ю.Г. Оптимизация и быстрое автоматическое дифференцирование. М.: ВЦ им. А.А. Дородницына РАН, 2013. 144 с.
Евтушенко Ю.Г., Зубов В.И. Об обобщенной методологии быстрого автоматического дифференцирования // Ж. вычисл. матем. и матем. физ. 2016. Т. 56. № 11. С. 1847–1862.
Албу А.Ф., Евтушенко Ю.Г., Зубов В.И. О выборе разностных схем при решении обратных коэффициентных задач // Ж. вычисл. матем. и матем. физ. 2020. Т. 60. № 10. С. 1643–1655.
Самарский А.А., Вабищевич П.Н. Вычислительная теплопередача. М.: Едиториал УРСС, 2003.
Douglas J., Rachford H.H. On the numerical solution of heat conduction problems in two and three space variables. Trans. Amer. Math. Soc. 1956. V. 82. P. 421–439.
Peaceman D.W., Rachford H.H. The Numerical Solution of Parabolic and Elliptic Differential Equations // Journal of the Society for Industrial and Applied Math. 1955. V. 3. № 1. P. 28–41.
Албу А.В., Зубов В.И. Опыт использования методологии Быстрого Автоматического Дифференцирования для решения обратных коэффициентных задач // Ж. вычисл. матем. и матем. физ. 2020. Т. 60. № 1. С. 18–28.
Самарский А.А. Теория разностных схем. М.: Наука, 1977.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал