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

Полный текст

ВВЕДЕНИЕ

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

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

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

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

В настоящей работе рассматривается одна из возможных постановок такой обратной коэффициентной задачи. Ее рассмотрение проводится на основе первой краевой задачи для одномерного нестационарного уравнения теплопроводности. Обратная коэффициентная задача сводится к вариационной задаче. В качестве целевого функционала выбрано среднеквадратичное отклонение рассчитанного теплового потока на границе образца от его экспериментального значения. Предложен алгоритм численного решения обратной коэффициентной задачи. В основе этого алгоритма лежит современная методология быстрого автоматического дифференцирования (БАД-методология, см. [15], [16]), позволившая успешно решить ряд сложных задач оптимального управления динамическими системами (см., например, [17]—[19]). Приведенные в работе примеры решения обратной коэффициентной задачи демонстрируют работоспособность предложенного алгоритма. В статье, кроме того, исследуется вопрос единственности получаемого решения сформулированной обратной задачи.

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

Рассматривается слой материала ширины L. Температура этого слоя в начальный момент времени известна. Также известно, как изменяется со временем температура на краях этого слоя. Распределение температурного поля в слое в каждый момент времени описывается решением следующей начально-краевой (смешанной) задачи:

$C\left(T\right)\frac{\partial T(x,t)}{\partial t}-\frac{\partial }{\partial x}\left(K\left(T\right)\frac{\partial T(x,t)}{\partial x}\right)=\mathrm{0,} (x,t)\in Q$, (1.1)

$T\left(x\mathrm{,0}\right)=w_0\left(x\right), (0\le x\le L)$, (1.2)

$T\left(\mathrm{0,}t\right)=w_1\left(t\right), T(L,t)=w_2\left(t\right), (0\le t\le \Theta )$. (1.3)

Здесь x — декартова координата точки в слое; t — время; $Q=\left\{(0<x<L)\times (0<t<\Theta )\right\}$; T(x, t) — температура вещества в точке с координатой x в момент времени t; C(T) — объемная теплоемкость вещества, K(T) — коэффициент теплопроводности; w₀(x) — заданная температура слоя в начальный момент времени; w₁(t) — заданная температура на левом краю слоя; w₂(t) — заданная температура на правом краю слоя.

Если зависимость объемной теплоемкости вещества C(T) и его коэффициента теплопроводности K(T) от температуры T известна, то, решая смешанную задачу (1.1) — (1.3), найдем распределение температуры T(x, t) в $\overline{Q}$. Задачу (1.1)—(1.3) ниже будем называть прямой задачей.

Если же зависимость объемной теплоемкости вещества и его коэффициента теплопроводности от температуры неизвестна, то представляет большой интерес задача определения этих зависимостей. Одной из возможных постановок такой задачи (эта задача относится к классу задач идентификации параметров модели) является следующая: найти такие зависимости объемной теплоемкости вещества C(T) и коэффициента теплопроводности K(T) от температуры T, при которых тепловые потоки ${\Pi }_0\left(t\right)=K\left(T\left(\mathrm{0,}t\right)\right)\cdot \frac{\partial T}{\partial x}\left(\mathrm{0,}t\right)$ на левом конце и ${\Pi }_L\left(t\right)=-K\left(T(L,t)\right)\cdot \frac{\partial T}{\partial x}(L,t)$ на правом конце, полученные в результате решения сформулированной задачи (1.1)-(1.3), мало отличаются от потоков P₀(t) и P_L(t) соответственно, полученных экспериментально. Мерой отклонения этих функций может служить величина (функционал «поток»):

$\Phi \left(C\left(T\right),K\left(T\right)\right)=\underset{0}{\overset{\Theta }{\int }}\left[{\beta }_0\left(t\right)\cdot {\left({\Pi }_0\left(t\right)-{{\rm P}}_0\left(t\right)\right)}^2+{\beta }_L\left(t\right)\cdot {\left({\Pi }_L\left(t\right)-{{\rm P}}_L\left(t\right)\right)}^2\right]dt$, (1.4)

где β₀(t) ≥ 0, β_L(t) ≥ 0 — заданные весовые параметры; P₀(t) и P_L(t) — заданные тепловые потоки на концах слоя. Таким образом, задача оптимального управления состоит в определении оптимального управления {C(T), K(T)} и соответствующего оптимального решения T(x, t) задачи (1.1)—(1.3), при котором функционал (1.4) достигает минимального значения.

Следует отметить, что задача идентификации объемной теплоемкости C(T) и коэффициента теплопроводности K(T), подобная описанной выше, рассматривалась в работе [20], только там в качестве целевого функционала выбирался функционал «поле»:

$\Phi \left(C\left(T\right),K\left(T\right)\right)=\underset{0}{\overset{\Theta }{\int }}\underset{0}{\overset{L}{\int }}{\left[T(x,t)-{\rm Y}(x,t)\right]}^2\cdot \mu (x,t)dxdt$, (1.5)

где μ(x, t) ≥ 0 — заданная весовая функция, T(x, t) — температурное поле, полученное в результате решения сформулированной задачи (1.1)—(1.3), Y(x, t) — температурное поле, полученное экспериментально. В [20] отмечено, что при использовании функционала «поле» решение задачи идентификации всегда неединственно и что целесообразно определять только отношение $\chi \left(T\right)=C\left(T\right)/K\left(T\right)$. В случае использования функционала «поток» (1.4) это утверждение оказывается несправедливым.

Получить аналитическое решение сформулированной обратной задачи (1.1)—(1.4) удается крайне редко. Поэтому необходимо разработать алгоритм численного решения сформулированной обратной задачи.

2. ЗАДАЧА ОПТИМАЛЬНОГО УПРАВЛЕНИЯ В ДИСКРЕТНОМ СЛУЧАЕ

При численном решении сформулированной задачи оптимального управления искомые функции C(T) и K(T), T ∈ [a, b] аппроксимировались непрерывными кусочно-линейными функциями следующим образом. Отрезок [a, b], на котором идентифицируются функции C(T) и K(T), разбивался точками ${\tilde{T}}_0=a, {\tilde{T}}_1, {\tilde{T}}_2, ..., {\tilde{T}}_N=b$ на N частей (они могут быть как равными, так и неравными). Каждой из точек ${\tilde{T}}_n(n=\mathrm{0,} ..., N)$ ставились в соответствие числа $c_n=C\left({\tilde{T}}_n\right), k_n=K\left({\tilde{T}}_n\right)$. Искомые функции C(T) и K(T) аппроксимировались непрерывными кусочно-линейными функциями с узлами в точках ${\left\{({\tilde{T}}_n, c_n)\right\}}_{n=0}^N$ и ${\left\{({\tilde{T}}_n, k_n)\right\}}_{n=0}^N$ так, что

$C\left(T\right)=c_n+\frac{c_{n+1}-c_n}{{\tilde{T}}_{n+1}-{\tilde{T}}_n}(T-{\tilde{T}}_n) при {\tilde{T}}_n\le T\le {\tilde{T}}_{n+1}, (n=\mathrm{0,...,}N-1)$, 

$K\left(T\right)=k_n+\frac{k_{n+1}-k_n}{{\tilde{T}}_{n+1}-{\tilde{T}}_n}(T-{\tilde{T}}_n) при {\tilde{T}}_n\le T\le {\tilde{T}}_{n+1}, (n=\mathrm{0,...,}N-1)$.

Область $Q=\left\{(x,t):(0<x<L)\times (0<t<\Theta )\right\}$ разбивалась сеточными линиями ${\left\{{\tilde{x}}_i\right\}}_{i=0}^I$ и ${\left\{{\tilde{t}}^j\right\}}_{j=0}^J$ на ряд прямоугольников с шагами $h_i={\tilde{x}}_{i+1}-{\tilde{x}}_i(i=\overline{\mathrm{0,}I-1}), {\tau }^j={\tilde{t}}^j-{\tilde{t}}^{j-1}(j=\overline{\mathrm{1,}J})$. В каждом узле $({\tilde{x}}_i,{\tilde{t}}^j)$ расчетной области Q, характеризуемом парой индексов (i, j), все функции задавались своими значениями в точке $({\tilde{x}}_i,{\tilde{t}}^j)$ (например, $T({\tilde{x}}_i,{\tilde{t}}^j)=T_i^j$). В каждом из полученных прямоугольников требовалось выполнение закона теплового баланса. В результате получалась хорошо известная конечно-разностная схема, аппроксимирующая смешанную задачу (1.1)—(1.3). Полученная система нелинейных алгебраических уравнений решалась итерационным методом с привлечением метода прогонки [21]. С помощью описанного подхода решалась смешанная задача (1.1)—(1.3) и находилась функция T(x, t) (точнее, ее аппроксимация $T_i^j$).

Целевой функционал (1.4) аппроксимировался функцией $F(c_0, c_1, ..., c_N, k_0, k_1, ..., k_N)$ конечного числа переменных следующим образом:

$\Phi \left(C\left(T\right),K\left(T\right)\right)\approx F=\sum _{j=1}^J\left[{\beta }_0\left(t^j\right)\cdot {\left({\Pi }_0\left(t^j\right)-{{\rm P}}_0\left(t^j\right)\right)}^2+{\beta }_L\left(t^j\right)\cdot {\left({\Pi }_L\left(t^j\right)-{{\rm P}}_L\left(t^j\right)\right)}^2\right]\cdot {\tau }^j$,

где

${\Pi }_0\left(t^j\right)=K\left(T_0^j\right)\cdot \left(T_1^j-T_0^j\right)/h_0, {\Pi }_L\left(t^j\right)=-K\left(T_I^j\right)\cdot \left(T_I^j-T_{I-1}^j\right)/h_{I-1}$.

Для минимизации целевого функционала в работе использовался градиентный метод. Хорошо известно, что для эффективной работы градиентных методов следует использовать точное значение градиента целевой функции. Вычислить значение градиента эффективно и с машинной точностью позволяет БАД-методология (см. [15], [16], [22]). Именно применение БАД-методологии позволило успешно решить ряд сложных задач оптимального управления. В работах [17]—[19], [22] использовался «ручной» вариант БАД-методологии, при котором все необходимые формулы выводились вручную и использовались в оптимизационном коде. Однако появились стандартные пакеты программ, в автоматическом режиме выполняющие примерно ту же работу (см. [23]—[25]). В статье [26] представлены результаты сравнения градиентов, вычисленных с помощью ручной версии БАД-методологии, с помощью стандартных пакетов программ, и сделан вывод, что стандартные пакеты программ позволяют практически с той же точностью вычислять градиенты сложных функций, что и ручная версия БАД-методологии. Поэтому в настоящей работе градиент целевой функции вычислялся с помощью пакета прикладного программного обеспечения Adept (см. [24]), при этом использовался обратный (Reverse, Adjoint) метод автоматического дифференцирования.

Для поиска минимума целевой функции применялся метод L-BFGS-B (см. [27]).

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

Минимальное значение целевого функционала (1.4), получаемое при решении сформулированной обратной задачи, зависит от заданных экспериментальных потоков P₀(t) и P_L(t). Если экспериментальные потоки принадлежат области достижимости функций Π₀(t) и Π_L(t), построенных в результате решения прямой задачи (1.1)—(1.3) при допустимых C(T) и K(T), то значение функционала равно нулю. В противном случае значение функционала больше нуля. При проведении численных экспериментов экспериментальные потоки P₀(t) и P_L(t) выбирались из области достижимости решений прямой задачи.

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

Все примеры первой группы базировались на том, что функция

$T_{*}(x,t)={\left(\frac{n+1}{\beta }\right)}^{\frac{1}{n+1}}\cdot {\left[\alpha \cdot (t+\gamma -x)\right]}^{\frac{1}{m-n}}, \alpha =\frac{(m-n)}{(n+1)}, \beta ={(n+1)}^{\frac{m-n}{m+1}}$

является аналитическим решением уравнения (1.1) при γ = const, C(T) = β · Tⁿ, K(T) = T^m.

В соответствии с этим использовались следующие входные данные задачи:

$$\begin{gathered} L=1, \gamma =1.5, Q=\left(\mathrm{0,1}\right)\times \left(\mathrm{0,}\Theta \right), \\ {w}_0\left(x\right)=T_{*}\left(x\mathrm{,0}\right)={\left(\frac{n+1}{\beta }\right)}^{\frac{1}{n+1}}\cdot {\left[\alpha \cdot (\gamma -x)\right]}^{\frac{1}{m-n}}, (0\le x\le 1), \\ {w}_2\left(t\right)=T_{*}\left(\mathrm{1,}t\right)={\left(\frac{n+1}{\beta }\right)}^{\frac{1}{n+1}}\cdot {\left[\alpha \cdot (t+\gamma -1)\right]}^{\frac{1}{m-n}}, (0\le t\le \Theta ), \\ {{\rm P}}_0\left(t\right)={\Pi }_0\left(t\right)=T_{*}^m\left(\mathrm{0,}t\right)\cdot \frac{\partial T_{*}}{\partial x}\left(\mathrm{0,}t\right), (0\le t\le \Theta ), \\ {{\rm P}}_L\left(t\right)={\Pi }_L\left(t\right)=-T_{*}^m(L,t)\cdot \frac{\partial T_{*}}{\partial x}(L,t), (0\le t\le \Theta ), \\ a=\underset{(x,t\in \overline{Q})}{\min }{T}_{*}(x,t), b=\underset{(x,t\in \overline{Q})}{\max }{T}_{*}(x,t). \end{gathered}$$ (3.1)

Что касается величины Θ, то она в расчетах принимала различные значения.

При решении прямой задачи использовалась равномерная сетка с параметрами I = 40 (количество интервалов вдоль оси x), J = 700 (количество интервалов вдоль оси t), которая обеспечивает достаточную точность вычисления поля температур и поля импульсов (сопряженных переменных). Отрезок [a, b] разбивался на 32 интервала (N = 32).

В первом примере первой группы расчетов рассматривалась задача нахождения объемной теплоемкости и коэффициента теплопроводности с входными данными (3.1) при m = 4 и n = 2. При этом полагалось, что в целевом функционале (1.4) весовые функции β₀(t), β_L(t) тождественно равны единице. Задача с такими параметрами исследовалась при разных значениях интервала времени Θ.

Вначале при решении задачи оптимизации функционала (1.4) интервал времени Θ = 1.5. Температура на левом конце слоя в этом случае изменяется от 1.246 до 1.762, а на правом — от 0.7192 до 1.438 (есть общий диапазон температур). В качестве начального приближения выбирались функции C(T) = 1.0 и K(T) = 1.0. Результаты полученного решения обратной задачи представлены на фиг. 1 (функция C(T)) и фиг. 2 (функция K(T)). Расчеты, выполненные с использованием других (отличных от констант) начальных приближений для функций C(T) и K(T), приводили к тому же решению обратной задачи.

Фиг. 1. Распределение объемной теплоемкости

Фиг. 2. Распределение коэффициента теплопроводности

Затем при решении той же задачи полагалось Θ = 0.5. Отметим, что здесь температура на левом конце слоя изменяется от 1.246 до 1.438, а на правом — от 0.7192 до 1.017 (есть диапазон температур, для которых потоки не задаются). В зависимости от начального приближения и параметров оптимизационного процесса получались различные решения. На фиг. 3—6 представлены два из получавшихся решений.

Фиг. 3. Распределение объемной теплоемкости

Фиг. 4. Распределение коэффициента теплопроводности

Фиг. 5. Распределение объемной теплоемкости

Фиг. 6. Распределение коэффициента теплопроводности

На фиг. 1—6 линии, помеченные как “Teor”, представляют собой ожидаемые решения C(T) = β · Tⁿ, K(T) = T^m обратной задачи (1.1)—(1.4) с данными (3.1), а линии, помеченные как “Opt”, — решения, полученные в результате минимизации целевого функционала (1.4).

Во втором примере первой группы расчетов рассматривалась задача нахождения объемной теплоемкости и коэффициента теплопроводности с входными данными (3.1) при m = 2 и n = 1. При этом, как и раньше, полагалось, что в целевом функционале (1.4) весовые функции β₀(t), β_L(t) тождественно равны единице. Задача с такими параметрами также исследовалась при разных значениях интервала времени Θ. Качественно ничего нового по сравнению с предыдущим примером здесь обнаружено не было: решение обратной задачи было единственным, если интервалы температур на левом и правом концах слоя (там, где задаются тепловые потоки) пересекаются, и решение неединственно, если на отрезке [a, b] имеется диапазон температур, для которых потоки не задаются. В предельном случае, когда интервалы температур на левом и правом концах слоя не пересекаются и непрерывно следуют один за другим (на левом конце слоя температура изменяется от 0.9449 до 1.5749, а на правом — от 0.3150 до 0.9449, Θ = 1.0) решение также оказывается единственным, но функционал стремится к нулевому значению чрезвычайно медленно. На фиг. 7 и фиг. 8 представлены промежуточные результаты решения обратной задачи, соответствующие значению целевого функционала (1.4) $\Phi \left(C\left(T\right),K\left(T\right)\right)\approx F\approx {10}^{-9}$.

Фиг. 7. Распределение объемной теплоемкости

Фиг. 8. Распределение коэффициента теплопроводности

Задача нахождения объемной теплоемкости и коэффициента теплопроводности с входными данными (3.1) решалась и при других значениях параметров m и n. На основе результатов, полученных при решении всех обратных задач первой группы, можно сделать вывод о том, что решение задачи определения объемной теплоемкости и коэффициента теплопроводности может быть неединственным.

В качестве первого примера второй группы рассматривалась обратная задача (1.1)—(1.4) со следующими входными данными:

$$\begin{gathered} L=\mathrm{1,} Q=\left(\mathrm{0,1}\right)\times \left(\mathrm{0,}\Theta \right), \\ {w}_0\left(x\right)=x^2, (0\le x\le 1), \\ {w}_1\left(t\right)=t, w_2\left(t\right)=t+\mathrm{1,} (0\le t\le \Theta ), \\ {{\rm P}}_0\left(t\right)={\Pi }_0\left(t\right)=\mathrm{0,} {{\rm P}}_L\left(t\right)={\Pi }_L\left(t\right)=-\mathrm{2,} (0\le t\le \Theta ), \\ a=\mathrm{0,} b=2, \end{gathered}$$ (3.2)

Заметим, что функция $T_{*}(x,t)=x^2+t$ является аналитическим решением уравнения (1.1) при C(T) = 2 и K(T) = 1, так что пара {C(T), K(T)} = {2, 1} и функция $T_{*}(x,t)=x^2+t$ являются оптимальным управлением и соответствующим оптимальным решением обратной задачи (1.1)—(1.4) с данными (3.2).

Величина Θ, как и ранее, принимала различные значения.

При времени протекания процесса Θ = 1.5 температура на левом конце слоя изменяется от 0.0 до 1.5, а на правом — от 1.0 до 2.5 (есть общий диапазон температур). В этом случае независимо от выбора начального приближения и параметров используемого оптимизационного алгоритма получалось одно и то же ожидаемое решение обратной задачи C(T) = 2, K(T) = 1.

Если же Θ = 0.5, то в этом случае температура на левом конце слоя изменяется от 0.0 до 0.5, а на правом — от 1.0 до 1.5 (есть диапазон температур, для которых потоки не задаются). Получить решение обратной задачи оказалось чрезвычайно сложно. Целевой функционал здесь имеет большое число локальных экстремумов, и используемые методы оптимизации первого порядка часто “застревали” в локальных минимумах. Удалось так подобрать параметры оптимизационного алгоритма, который приводил к нахождению ожидаемого решения обратной задачи. При этом следует отметить, что элементы сходящейся к решению последовательности имели довольно экзотичный вид (см. фиг. 9, фиг. 10).

Фиг. 9. Распределение объемной теплоемкости

Фиг. 10. Распределение коэффициента теплопроводности

Во втором примере второй группы рассматривалась обратная задача (1.1)—(1.4) с такими входными данными:

$$\begin{gathered} L=\mathrm{1,} \Theta =\mathrm{1,} Q=\left(\mathrm{0,1}\right)\times \left(\mathrm{0,1}\right), \\ {w}_0\left(x\right)=9/(x+1), (0\le x\le 1), \\ {w}_1\left(t\right)=9/(5t+1), (0\le t\le 1), \\ {w}_2\left(t\right)=9/(7t+2), (0\le t\le 1), \\ {{\rm P}}_0\left(t\right)={\Pi }_0\left(t\right)=K\left(T_{**}\left(\mathrm{0,}t\right)\right)\cdot \frac{\partial T_{**}}{\partial x}\left(\mathrm{0,}t\right), (0\le t\le \Theta ), \\ {{\rm P}}_L\left(t\right)={\Pi }_L\left(t\right)=-K\left(T_{**}(L,t)\right)\cdot \frac{\partial T_{**}}{\partial x}(L,t), (0\le t\le \Theta ), \\ a=\mathrm{1,} b=9, \end{gathered}$$

где экспериментальное поле T_**(x, t) строилось численно как решение следующей смешанной задачи (1.1)—(1.3):

$$\begin{gathered} C\left(T\right)\frac{\partial T(x,t)}{\partial t}-\frac{\partial }{\partial x}\left(K\left(T\right)\frac{\partial T(x,t)}{\partial x}\right)=\mathrm{0,} (0<x<1), (0<t<1), \\ T\left(x\mathrm{,0}\right)=9/(x+1), (0<x<1), \\ T\left(\mathrm{0,}t\right)=9/(5t+1), (0<t<1), \\ T(L,t)=9/(7t+2), (0<t<1), \end{gathered}$$

при

$$\begin{gathered} K\left(T\right)=\left\{\begin{array}{ll}0.1\cdot (T-3)\cdot (T-6)\cdot (T-7)+3.4,& T\ge 3,\\ 1.2\cdot (T-3)+3.4,& T<3,\end{array}\right. \\ C\left(T\right)=4\cdot T^2, T\in [a,b], a=\mathrm{1,} b=9. \end{gathered}$$

Отметим, что в этом примере отрезки изменения температуры на левом и правом краях слоя (там, где задаются тепловые потоки) имеют непустое пересечение.

Здесь, как и в первом примере второй группы, обратная задача решалась численно при разных начальных приближениях. Независимо от начального приближения получалось одно и то же решение. На фиг. 11 и фиг. 12 представлены результаты решения задачи при начальном приближении C(T) = 4, K(T) = 1.

Фиг. 11. Распределение объемной теплоемкости

Фиг. 12. Распределение коэффициента теплопроводности

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

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

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

Искомые объемная теплоемкость C(T) и коэффициент теплопроводности K(T) должны принадлежать области определения, т. е. удовлетворять условиям $C\left(T\right)\in C\left(\right[a,b\left]\right)$, $K\left(T\right)\in C^1\left(\right[a,b\left]\right)$, C(T) > 0, K(T) > 0 для всех T ∈ [a, b]. Под областью достижимости будем понимать множество функций $T(x,t)\in C_{x,t}^{\mathrm{2,1}}\left(Q\right)\cap C_{x,t}^{\mathrm{1,0}}\left(\overline{Q}\right)$, ${{\rm P}}_0\left(t\right)\in C\left(\left[\mathrm{0,}\Theta \right]\right)$, ${{\rm P}}_L\left(t\right)\in C\left(\left[\mathrm{0,}\Theta \right]\right)$ таких, что функции T(x, t) удовлетворяют заданным условиям (1.2),(1.3), для каждой из функций T(x, t) в области определения найдутся функции C(T) и K(T), при которых T(x, t) является решением смешанной задачи (1.1)—(1.3), а функции P₀(t) и P_L(t) определяются соотношениями

${{\rm P}}_0\left(t\right)=K\left(T\left(\mathrm{0,}t\right)\right)\cdot \frac{\partial T}{\partial x}\left(\mathrm{0,}t\right), {{\rm P}}_L\left(t\right)=-K\left(T(L,t)\right)\cdot \frac{\partial T}{\partial x}(L,t)$.

Будем под τ-изолинией функции T(x, t) будем понимать, как это принято, множество точек (x, t) области Q, в которых значение функции T(x, t) имеет одно и то же значение τ. При этом отметим, что изолиния (множество) может представлять собой объединение нескольких, не связанных между собой линий.

Теорема. «Пусть функции T(x, t), P₀(t) и P_L(t) принадлежат области достижимости, $a=\underset{(x,t\in \overline{Q})}{\min }T(x,t)$, $b=\underset{(x,t\in \overline{Q})}{\max }T(x,t)$. Пусть на отрезке [a, b] имеется интервал (γ, δ) такой, что $w_1\left(t\right)=T\left(\mathrm{0,}t\right)\notin (\gamma ,\delta )$, $w_2\left(t\right)=T(L,t)\notin (\gamma ,\delta )$ для всех 0 ≤ t ≤ Θ. Если найдется отрезок $[\alpha ,\beta ]\in (\gamma ,\delta )\in [a,b]$ такой, что на каждой τ-изолинии отрезка τ ∈ [α, β] функции $\frac{\partial T(x,t)}{\partial t}$ и ${\left(\frac{\partial T(x,t)}{\partial x}\right)}^2$, рассматриваемые как функции некоторого параметра s вдоль этой изолинии, являются линейно зависимыми и функция $\frac{\partial T(x,t)}{\partial t}$ не равна тождественно нулю, то решение {C(T), K(T)} обратной задачи (1.1)—(1.4) неединственно)».

Доказательство. На первом этапе доказательства теоремы обратим внимание на линейную зависимость функций $\frac{\partial T(x,t)}{\partial t}$ и ${\left(\frac{\partial T(x,t)}{\partial x}\right)}^2$ вдоль каждой τ-изолинии отрезка τ ∈ [α, β]. Это означает, что существуют функции A(T) и B(T), T ∈ [α, β], такие что

$A\left(T\right)\cdot \frac{\partial T(x,t)}{\partial t}+B\left(T\right)\cdot {\left(\frac{\partial T(x,t)}{\partial x}\right)}^2=0$,

причем $A^2\left(T\right)+B^2\left(T\right)>0$. Более того, A(T) ≠ 0 для всех T ∈ [α, β]. Действительно, предположив, что А(τ) = 0 для какого-то значения τ ∈ [α, β], придем к противоречию с условием отличия от тождественного нуля производной $\frac{\partial T(x,t)}{\partial t}$. Таким образом, справедливо следующее равенство:

$\frac{\partial T(x,t)}{\partial t}=\tilde{B}\left(T\right)\cdot {\left(\frac{\partial T(x,t)}{\partial x}\right)}^2$, (4.1)

где

$\tilde{B}\left(T\right)=-\frac{B\left(T\right)}{A\left(T\right)}\in C\left(\right[\alpha ,\beta \left]\right)$.

Так как функции T(x, t), P₀(t) и P_L(t) по предположению принадлежат области достижимости, то существуют функции C_*(T) и K_*(T) из области определения, при которых функция T(x, t) является решением смешанной задачи (1.1)—(1.3), а функции P₀(t) и P_L(t) определяются равенствами

${{\rm P}}_0\left(t\right)=K_{*}\left(T\left(\mathrm{0,}t\right)\right)\cdot \frac{\partial T}{\partial x}\left(\mathrm{0,}t\right), {{\rm P}}_L\left(t\right)=-K_{*}\left(T(L,t)\right)\cdot \frac{\partial T}{\partial x}(L,t)$.

Покажем, что в области определения кроме них найдутся другие функции C(T) и K(T), при которых эта же функция T(x, t) является решением смешанной задачи (1.1)—(1.3) и

${{\rm P}}_0\left(t\right)=K\left(T\left(\mathrm{0,}t\right)\right)\cdot \frac{\partial T}{\partial x}\left(\mathrm{0,}t\right), {{\rm P}}_L\left(t\right)=-K\left(T(L,t)\right)\cdot \frac{\partial T}{\partial x}(L,t)$,

Пусть $\psi \left(T\right)\in \stackrel{0}{C}\left(\right[\alpha ,\beta \left]\right)$ — произвольная функция такая, что $\underset{\alpha }{\overset{\beta }{\int }}\tilde{B}\left(\xi \right)\cdot \psi \left(\xi \right)d\xi =0$. Определим функции $f\left(T\right)\in C\left(\right[\alpha ,\beta \left]\right)$ и $\sigma \left(T\right)\in C^1\left(\right[\alpha ,\beta \left]\right)$ соотношениями

$f\left(T\right)=\frac{{K^{\text{'}}}_{*}\left(T\right)}{K_{*}\left(T\right)}+\tilde{B}\left(T\right)\cdot \psi \left(T\right), \sigma \left(T\right)=\exp \left[\underset{\alpha }{\overset{T}{\int }}f\left(\xi \right)d\xi \right]$.

Построим новые объемную теплоемкость C(T) и коэффициент теплопроводности K(T)

$K\left(T\right)=\left\{\begin{array}{ll}{K}_{*}\left(T\right),& a\le T<\alpha ,\\ K_{*}\left(\alpha \right)\cdot \sigma \left(T\right),& \alpha \le T\le \beta ,\\ K_{*}\left(T\right)\cdot K_{*}\left(\alpha \right)\cdot \sigma \left(\beta \right)/K_{*}\left(\beta \right),& \beta <T\le b,\end{array}\right.$ (4.2)

$C\left(T\right)=\left\{\begin{array}{ll}{C}_{*}\left(T\right),& a\le T<\alpha ,\\ K\left(T\right)\cdot \left(C_{*}\left(T\right)/K_{*}\left(T\right)+\psi \left(T\right)\right),& \alpha \le T\le \beta ,\\ C_{*}\left(T\right)\cdot K\left(T\right)/K_{*}\left(T\right),& \beta <T\le b,\end{array}\right.$

и покажем, что они также являются решением обратной задачи для рассматриваемых функции-поля T(x, t) и потоков P₀(t), P_L(t) из области достижимости.

Во-первых, непосредственной проверкой нетрудно убедиться в том, что построенные функции $K\left(T\right)\in C^1\left(\right[a,b\left]\right)$, $C\left(T\right)\in C\left(\right[a,b\left]\right)$ и что при всех T ∈ [α, β] справедливо соотношение

$\frac{K^{\text{'}}\left(T\right)}{K\left(T\right)}-\tilde{B}\left(T\right)\cdot \frac{C\left(T\right)}{K\left(T\right)}=\frac{{K^{\text{'}}}_{*}\left(T\right)}{K_{*}\left(T\right)}-\tilde{B}\left(T\right)\cdot \frac{C_{*}\left(T\right)}{K_{*}\left(T\right)}$, (4.3)

а при всех T ∈ [a, α] и T ∈ [β, b] — соотношения

$\frac{K^{\text{'}}\left(T\right)}{K\left(T\right)}=\frac{{K^{\text{'}}}_{*}\left(T\right)}{K_{*}\left(T\right)}, \frac{C\left(T\right)}{K\left(T\right)}=\frac{C_{*}\left(T\right)}{K_{*}\left(T\right)}$. (4.4)

Во-вторых, используя соотношения (4.3), (4.4), покажем, что при заданной в условии теоремы функции T(x, t) уравнение (1.1) выполняется не только при заданных (C_*(T), K_*(T)), но и при вновь построенных (C(T), K(T)). Действительно, при предположениях, сформулированных в условиях теоремы, уравнение (1.1) можно представить в виде

$C\left(T\right)\cdot \frac{\partial T(x,t)}{\partial t}-K^{\text{'}}\left(T\right)\cdot {\left(\frac{\partial T(x,t)}{\partial x}\right)}^2-K\left(T\right)\cdot \frac{{\partial }^{2}T(x,t)}{\partial x^2}=0$. (4.5)

В тех точках области Q, в которых T ∈ [a, α] или T ∈ [β, b], имеем

$\begin{array}{c}\Psi (x,t)\equiv C\left(T\right)\cdot \frac{\partial T(x,t)}{\partial t}-K^{\text{'}}\left(T\right)\cdot {\left(\frac{\partial T(x,t)}{\partial x}\right)}^2-K\left(T\right)\cdot \frac{{\partial }^{2}T(x,t)}{\partial x^2}=\\ =K\left(T\right)\cdot \left[\frac{C\left(T\right)}{K\left(T\right)}\cdot \frac{\partial T(x,t)}{\partial t}-\frac{K^{\text{'}}\left(T\right)}{K\left(T\right)}\cdot {\left(\frac{\partial T(x,t)}{\partial x}\right)}^2-\frac{{\partial }^{2}T(x,t)}{\partial x^2}\right]=\\ =K\left(T\right)\cdot \left[\frac{C_{*}\left(T\right)}{K_{*}\left(T\right)}\cdot \frac{\partial T(x,t)}{\partial t}-\frac{{K^{\text{'}}}_{*}\left(T\right)}{K_{*}\left(T\right)}\cdot {\left(\frac{\partial T(x,t)}{\partial x}\right)}^2-\frac{{\partial }^{2}T(x,t)}{\partial x^2}\right]=\\ =\frac{K\left(T\right)}{K_{*}\left(T\right)}\cdot \left[C_{*}\left(T\right)\cdot \frac{\partial T(x,t)}{\partial t}-{K^{\text{'}}}_{*}\left(T\right)\cdot {\left(\frac{\partial T(x,t)}{\partial x}\right)}^2-K_{*}\left(T\right)\cdot \frac{{\partial }^{2}T(x,t)}{\partial x^2}\right]=0.\end{array}$

Здесь использовались соотношения (4.4).

В тех же точках области Q, в которых T ∈ [α, β], имеем

$\begin{array}{c}\Psi (x,t)\equiv C\left(T\right)\cdot \frac{\partial T(x,t)}{\partial t}-K^{\text{'}}\left(T\right)\cdot {\left(\frac{\partial T(x,t)}{\partial x}\right)}^2-K\left(T\right)\cdot \frac{{\partial }^{2}T(x,t)}{\partial x^2}=\\ =K\left(T\right)\cdot \left[\frac{C\left(T\right)}{K\left(T\right)}\cdot \frac{\partial T(x,t)}{\partial t}-\frac{K^{\text{'}}\left(T\right)}{K\left(T\right)}\cdot {\left(\frac{\partial T(x,t)}{\partial x}\right)}^2-\frac{{\partial }^{2}T(x,t)}{\partial x^2}\right]=\end{array}$ (используем условие (4.3))

$=K\left(T\right)\cdot \left\{\frac{C\left(T\right)}{K\left(T\right)}\cdot \frac{\partial T(x,t)}{\partial t}-\left[\frac{{K^{\text{'}}}_{*}\left(T\right)}{K_{*}\left(T\right)}-\tilde{B}\left(T\right)\cdot \left(\frac{C_{*}\left(T\right)}{K_{*}\left(T\right)}-\frac{C\left(T\right)}{K\left(T\right)}\right)\right]\cdot {\left(\frac{\partial T(x,t)}{\partial x}\right)}^2-\frac{{\partial }^{2}T(x,t)}{\partial x^2}\right\}=$ (используем условие (4.1) в первой квадратной скобке)

$=K\left(T\right)\cdot \left\{\frac{C_{*}\left(T\right)}{K_{*}\left(T\right)}\cdot \tilde{B}\left(T\right)\cdot {\left(\frac{\partial T(x,t)}{\partial x}\right)}^2-\frac{{K^{\text{'}}}_{*}\left(T\right)}{K_{*}\left(T\right)}\cdot {\left(\frac{\partial T(x,t)}{\partial x}\right)}^2-\frac{{\partial }^{2}T(x,t)}{\partial x^2}\right\}=$ (используем условие (4.1) в первом слагаемом)

$\begin{array}{c}=K\left(T\right)\cdot \left\{\frac{C_{*}\left(T\right)}{K_{*}\left(T\right)}\cdot \frac{\partial T(x,t)}{\partial t}-\frac{{K^{\text{'}}}_{*}\left(T\right)}{K_{*}\left(T\right)}\cdot {\left(\frac{\partial T(x,t)}{\partial x}\right)}^2-\frac{{\partial }^{2}T(x,t)}{\partial x^2}\right\}=\\ =\frac{K\left(T\right)}{K_{*}\left(T\right)}\cdot \left[C_{*}\left(T\right)\cdot \frac{\partial T(x,t)}{\partial t}-{K^{\text{'}}}_{*}\left(T\right)\cdot {\left(\frac{\partial T(x,t)}{\partial x}\right)}^2-K_{*}\left(T\right)\cdot \frac{{\partial }^{2}T(x,t)}{\partial x^2}\right]=0.\end{array}$

Таким образом, прямая задача (1.1)—(1.3), имеет одно и то же решение T(x, t) как для функций C_*(T) и K_*(T), так и для семейства функций C(T) и K(T).

В-третьих, покажем, что при всех T ∈ [a, α] и T ∈ [β, b] функции K(T) совпадают с функцией K_*(T). Что касается отрезка T ∈ [a, α], то здесь совпадение следует из определения функций K(T) (см. (4.2)). Для T ∈ [β, b] имеем

$\begin{array}{c}\sigma \left(\beta \right)=\exp \left[\underset{\alpha }{\overset{\beta }{\int }}\left(\frac{{K^{\text{'}}}_{*}\left(\xi \right)}{K_{*}\left(\xi \right)}+\tilde{B}\left(\xi \right)\cdot \psi \left(\xi \right)\right)d\xi \right]=\exp \left[\underset{\alpha }{\overset{\beta }{\int }}\left(\frac{{K^{\text{'}}}_{*}\left(\xi \right)}{K_{*}\left(\xi \right)}\right)d\xi \right]=\\ =\exp \left[\underset{\alpha }{\overset{\beta }{\int }}d\ln \left(K_{*}\left(\xi \right)\right)\right]=\exp \left[\ln \frac{K_{*}\left(\beta \right)}{K_{*}\left(\alpha \right)}\right]=\frac{K_{*}\left(\beta \right)}{K_{*}\left(\alpha \right)}.\end{array}$

Отсюда получим

$K\left(T\right)=K_{*}\left(T\right)\cdot K_{*}\left(\alpha \right)\cdot \sigma \left(\beta \right)/K_{*}\left(\beta \right)=K_{*}\left(T\right)\cdot K_{*}\left(\alpha \right)\cdot \frac{K_{*}\left(\beta \right)}{K_{*}\left(\alpha \right)}/K_{*}\left(\beta \right)=K_{*}\left(T\right)$.

Как следует из соотношений (4.2), при T ∈ [a, α] и T ∈ [β, b] функции C(T) также совпадают с функцией C_*(T).

Итак, при выполнении условий, сформулированных в теореме, обратная задача (1.1)—(1.4) наряду с решением C_*(T) и K_*(T) имеет семейство решений C(T) и K(T), определяемое формулами (4.2).

Замечание. Примерами, иллюстрирующими теорему, могут служить примеры из первой группы расчетов при Θ = 0.5.

ЗАКЛЮЧЕНИЕ

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

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

Об авторах

А. Ю. Горчаков

Email: vladimir.zubov@mail.ru
Россия, Москва

В. И. Зубов

Автор, ответственный за переписку.
Email: vladimir.zubov@mail.ru
Россия, Москва

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

Дополнительные файлы

Доп. файлы

Действие

1. JATS XML

Скачать

2. Фиг. 1. Распределение объемной теплоемкости

Скачать (90KB)

Метаданные

3. Фиг. 2. Распределение коэффициента теплопроводности

Скачать (93KB)

Метаданные

4. Фиг. 3. Распределение объемной теплоемкости

Скачать (82KB)

Метаданные

5. Фиг. 4. Распределение коэффициента теплопроводности

Скачать (85KB)

Метаданные

6. Фиг. 5. Распределение объемной теплоемкости

Скачать (93KB)

Метаданные

7. Фиг. 6. Распределение коэффициента теплопроводности

Скачать (92KB)

Метаданные

8. Фиг. 7. Распределение объемной теплоемкости

Скачать (69KB)

Метаданные

9. Фиг. 8. Распределение коэффициента теплопроводности

Скачать (71KB)

Метаданные

10. Фиг. 9. Распределение объемной теплоемкости

Скачать (104KB)

Метаданные

11. Фиг. 10. Распределение коэффициента теплопроводности

Скачать (118KB)

Метаданные

12. Фиг. 11. Распределение объемной теплоемкости

Скачать (88KB)

Метаданные

13. Фиг. 12. Распределение коэффициента теплопроводности

Скачать (88KB)

Метаданные