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

Вычислительные технологии

Том 18, № 1, 2013

Метод численного решения обратных нелинейных

задач по восстановлению компонентов тензора

>к

теплопроводности анизотропных материалов*

С. А. Колесник Московский авиационный институт, (Национальный исследовательский университет), Россия e-mail: [email protected]

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

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

Введение

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

— нестационарное температурное поле многомерно по пространственным переменным;

— вектор плотности теплового потока не ортогонален изотермам, вследствие чего все координатные направления равнозначны и главное из них выделить невозможно;

— дифференциальное уравнение теплопереноса содержит смешанные производные, что не позволяет разделить переменные по координатным направлениям;

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

* Работа выполнена при финансовой поддержке ФЦП "Научные и научно-педагогические кадры инновационной России" на 2009-2013 годы, госконтракт № 16.740.11.0454 и грантов РФФИ 12-01-31231 мол_а, 12-01-33095 мол_а_вед.

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

Обратные задачи теплопереноса для изотропных материалов рассматривались ранее в работах [1-3]. В [4] решена обратная задача теплопроводности в анизотропном полупространстве на основе аналитического решения, полученного в [5], в [6] решались задачи по восстановлению источника для уравнения диффузии.

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

1. Постановка задачи

В прямоугольной пластине из анизотропного материала 1\ х /2 (рис. 1) рассматривается следующая коэффициентная обратная задача теплопереноса по определению компонентов тензора теплопроводности Лц (Т) , Л12 (Т) , Л21 (Т) , Л22 (Т):

дХ (Ли (Т) §) + | (Л12 (Т) %) + | (Л.21 (Т) ) + | (Л22 (Т) |) =

Х е (0; /1), у е (0; /2), ¿> 0, (1)

Т(х, 0,£) = Т(Х,М) = Т(0,у,¿) = Т(/1 ,у,г) = Т|г , (2)

Т(х,у, 0) = То. (3)

Рис. 1. Расчётная область; Ох, Оу — оси декартовой системы координат, О£, Оп — главные оси тензора теплопроводности, ориентированные относительно декартовой системы координат под углом р

При этом граничное условие принимается в качестве максимального значения температур T |г = Tmax, а начальное — в качестве минимального T0 = Tmin, т. е.

T ■ < T < T (4)

± min < ± < ^max' V 7

Для замыкания коэффициентной обратной задачи теплопереноса в анизотропной пластине задаются экспериментальные значения температур в девяти точках (см. рис. 1) в зависимости от времени

T ((x,y)i ) = Tifc, i = 1,9, k = IK (5)

Нелинейные компоненты тензора теплопроводности определяются через главные коэффициенты Ag, An и угол ориентации ^ главных осей и On следующим образом [5]:

Ац (T) = Ag (T) cos2 + An (T) sin2 A22 (T) = Ag (T) sin2 + An (T) cos2 A12 (T) = A21 (T) = (Ag (T) - An (T)) sin ^ cos ^ (6)

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

Одним из эффективных методов решения нелинейных коэффициентных обратных задач для уравнений параболического типа вообще и анизотропного переноса тепла, в частности, является метод параметрической идентификации, [3, 4], в котором искомые функции A11 (T) , A12 (T) = A21 (T) , A22 (T) находятся в виде линейной комбинации базисных функций Nm (T) , задаваемых на конечных элементах — конечных отрезках ATm, m = I, M — I (Tmin < T < Tmax) , причём эти базисные функции приписаны

m _

каждому узлу Tm = Tmin + ATl-1 (AT0 = 0), m = I, M, и ортогональны на отрезке

i=1

T E [Tmin, Tmax] в пространстве L2. Здесь используются линейно непрерывные базисные функции [3, 7]

Nm (T) =

Нелинейные компоненты тензора теплопроводности, зависящие от температуры, в методе параметрической идентификации представляются в виде линейных комбинаций базисных функций (Т)

м

Л11 (Т) « ^ А^ (Т), (8)

т=1

T — Tm-1

T — T :

T m T m—1 Tm+1 — T

T — T :

Tm+1 Tm

T < Tm—1, Tm—1 < T < Tm,

Tm < T < Tm+1,

T > T

max

m = I,M.

0

м

Л22 (Т) « Е (Т), (9)

т=1

м

Л12 (Т) « Е Л^ (Т) , (10)

т=1

где коэффициенты линейных комбинаций Л™, Лт2 на каждом т-м конечном элементе т =1, (М — 1) являются искомыми величинами.

На основании принципа максимума можно утверждать, что температура внутри области будет удовлетворять неравенству Тт;п < Т(ж,у, ¿) < Ттах. Для дальнейших рассуждений введём следующие обозначения: Л^ = Л11 (Тт;п), Л^ = Л12 (Тт;п), Л^2 =

Л22 (Ттт) , . . . , Л11 = Л11 (Ттах) , Л12 = Л12 (Ттах)) Л22 = Л22 (Ттах).

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

Для определения постоянных компонентов вектора Л = (Л;[1,..., Лм1, Л1!2,..., Л12!, Л12,..., Лм)т, т =1, М, в выражениях (8)-(10) вводится квадратичный функционал

1 9 К о

5 (Л) = 2 ЕЕ К* (Л) —

2

г=1 *=1

в виде суммы по пространственно-временным переменным квадратов отклонения экспериментальных значений Т^ в точках ((ж, у^, ¿*) от расчётных Т^* (Л) = ((ж, у)^ Ьк, Л. В случае отсутствия экспериментальных значений температур в качестве последних принимаются результаты численного решения по приемлемым характеристикам Л11 (Т) , Л22 (Т) , Л12 (Т) , считающимся искомыми. При этом в экспериментальные значения может быть добавлена относительная 8 либо абсолютная А погрешность. Предполагается, что при достижении стационарного значения функционала (11) искомые характеристики, заложенные в экспериментальных значениях Т^*, приближённо совпадут с характеристиками, по которым получены расчётные значения температур.

Для минимизации функционала используется неявный метод градиентного спуска

Л(п+1) = Л(п) — а^гаа 5 (Л(п+1)) , (12)

где п — номер итерации, ап — параметрические шаги, выбираемые достаточно малыми с подчинением условию (ап > 0)

5 (Л(п+1)) < 5 (Л(п)). (13)

По условию (13) первоначальное значение а0 может быть выбрано произвольно, например а0 = 0.01. Тогда, если в результате следующей итерации условие (13) не выполнилось, то ап на этой итерации уменьшается и расчёт на ней повторяется, в противном случае — при выполнении (13) для следующей итерации — ап увеличивается. Окончание итерационного процесса устанавливается по близости к нулю grad 5 (Л(п+1)) , т.е. при выполнении условия

^ 5 (Л(п+1))|< е, (14)

где е — заданная точность.

2

Для вычисления градиента функционала в (11) с последующей подстановкой его компонентов в (12) и определения вектора АА(п) = А(га+1) — А(п) разложим в ряд Тейлора функцию Тг,& (А(га+1)) в окрестности А(п), сохраняя линейные относительно АА(п) члены. В результате получим

5 (А(

га+1)

9 Ко

Е Е

г=1 к=1

3М

Т,* А(п)

1=1

дТ^ (А( дА/

га)

АА(п) |— Та

+ 0(||АА|Г

'15)

Компоненты градиента функционала (15) имеют вид

д5 (А(п+1)) _Е Е

дА/

х

г=1 к=1

дТ^ (А(п) дА/

(Ч* (А(п)) — Т,*) + ^ дТ-,;АА(га^ АА((п)

к=1

X

+

А

дА

3м

£ дТ^(А(П)) АА(га)

9 Ко

Е Е

г=1 *=1

./=1 3м

дА/

Т^ (А(га)) — ТгИ

дТ^ (А(п))

/=1

дА/

АА

(га)

дТгЛ (А(п))

дА/

I = 1, 3М,

'16)

где 3М — количество неизвестных параметров.

Представим (16) в следующей векторно-матричной форме:

grad 5 (А(п+1)) = (А(п)) (т (А(п)) — т) + (А(п)) £ (А(п)) АА(п), (17

£ (А(

га)

( и1 ((ж,у)1,*1,А(га)) ... V1 ((ж,у)1,*1,А(га)) ... ^1((ж,у)1,¿1, А(п)) ... ^

и1 ((ж,у)2,ЛА(га)) ... V1 ((ж,у)2,^1,А(га)) ... ш1 ((ж,у)2, ¿1, А(п)) ...

18)

V«1 ((х,у)9, А(п)) ... V1 ((ж,у)9,£Ко, А(п)) ... ^ ((х,у)9,^Ко, А(п)) ... /

Элементы матрицы (18) (её можно назвать матрицей коэффициентов чувствительности) в точках ((ж,у)^, определяются из решения сопряжённых задач относительно производных от прямой задачи (1)-(4) по каждому компоненту А/, I = 1,.., 3М, и имеют смысл коэффициентов чувствительности температуры при изменении параметров А/:

ит((х,уМк, А)

дТ ((ж,уМ* ,А) дАТ1

, А)

Векторы Т (А(га)) — Т и АА(п) в (17) имеют вид

, , А) =

дТ ((х,уМ*,А)

дТ((х,у)г,^*, А

дАт ,

дАт2

'19)

Т (А(п)) — Т = (Тм (А(п)) — Тм) ,..., (Т/,1 (А(п)) — Т/д) , (Т1,2 (А(п)) — Т^) ,..., (Т/,2 (А(п)) — Т/,2) ,..., (Т/,Ко (А(п)) — Т^)

т

2

1

2

АЛ(п) = АЛ

йп),

.., АЛМ1(п), АЛ12п),..., АЛ22(п), АЛ12п),..., АЛМ2(п)

т

(21)

Подставляя (17) в (12), получим

АЛ(п) = —ап (Л(п)) (Т (Л(п)) — т) + (Л(п)) £ (Л(п)) АЛ(п)

откуда

(Е + ап^т (Л(п)) £ (Л(п))) АЛ(п) = —ап^т (Л(п)) (т (Л(п)) — т)

или

АЛ(п) = —ап (Е + ап^т (Л(п)) £ (Л(п)))-1 (Л(п)) (т (Л(п)) — т)

(22)

Для определения расчётных значений температур Т^* (Л) , входящих в функционал (11), используется экономичный абсолютно устойчивый метод переменных направлений с экстраполяцией (МПНЭ), подробно изложенный и обоснованный в [8, 9], а для нахождения элементов матрицы коэффициентов чувствительности необходимо решить 3М-независимые начально-краевые задачи относительно функций

и1(ж, у, *)

дТ

злЦ !

им (ж, у, *)

дТ

длм;

^1(ж, у, *)

дТ дТ

дл|2;

.Vм (ж, у,*) =

, (ж, у, *) =

дТ

длм,

дТ

дЛм,

которые можно получить, продифференцировав задачу (1)-(4) по соответствующим параметрам Л1ь...,Лм, Л!2,..., Л^^, Л12,...,Лм.

Продифференцируя уравнения задачи (1)-(4) по параметру Л^, получим

д_ дж

ду

м

м

луг ) + £ ЛТ1 дж + (£ ЛГ^„(Г)

т=1

11

ут=1

д_ дж

А

ду

м

Е л:

м

, дЛт(ТЛ 5Т

12 дЛ?1 ) ду

м

+ £ лт,лт(Т)

Л

Е Л52£ + (Е лВл„(Г)

ут=1

м

д дТ дж дЛС1

д дТ

+

дЖт (Т) \ дТ ду

дж м

ду дЛи д дТ

дждлЦ

дЛ?1

Л ЕГ Л^Т) ) (

д дТ

V ду длЦ

ут=1 / 4 ^ 11

ж е [0;/1] , у е [0; /2] , *> 0

= ср

+ +

д ( дТ

дА дЛ11

дТ . дТ дТ . . дТ

(ж, М) = тггр- (ж,/2,*) = тггр- (0,у,^) = тггр- (/1,у,*) = °

дЛ?1

дЛ?1

дЛ?1

дЛп

дТ

дЛЦ

(ж, у, 0) = 0.

(23)

(24)

(25)

дТ ~ м

Произведём замену в (23)-(25): плр = ир(х,у,*), А11 (Т) = ^ А^Ж^,(Т). В резуль-

т=1

дА?1

тате получим следующую начально-краевую дифференциальную задачу относительно неизвестной функции ир(х,у,*):

д_ дх

м

жр(Т ) + V Ат

дЖт(ТЛ дТ

т=1

дА1

11

дх + Л11 (Т) V дх

д«р

+

д_ дх

А

дУ

м

ЕаГ2дТ + Л.2 (Т/д«р

ут=1

м

дА?1

дУ

дУ

ЕАЗЗД дХ + Л, (Т)(

+

+

д_ ду

м

А

ут=1

д^т(ТЛ дТ + Л (Т^д«р

дАр1 ) + Л22 (ТЧ

= ср-

д«р

д*

х е [0; /1], у е [0; /2], 0, ир(х, 0,*) = ир(х,/2,*) = «р(0, у, *) = ир(/1,у,*) = 0, ир(х,у, 0) = 0.

(26) (27)

Аналогично продифференцировав задачу (1)-(4) по параметрам А^2, А222 и введя обозначения дТ/дАр2 = (х, у,*) и дТ/дА^2 = (х, у,*), получим начально-краевые дифференциальные задачи относительно неизвестных функций ^р(х,у,*) и ^р(х,у,*) соответственно:

А

дх

Ат д^т (Т) \ дТ + Л (т/ д^р а11^Тр I + А11 (Т )

ут=1

дА?2

дх

дх

+

А

дх

А

дУ

м

N (Т) + £ А

дЖт(ТЛ дТ

т_

12 дАр2

^ + Л12 (Т) —

дУ

-

дУ /

т=1

N (т) + Е А?2 ^) % + Л, (Т) (Щ)

+

+

А

дх

м

дЖт(Т) \ дТ + Л (Т/д^р А22 - I ~ + А22 (Т )

ут=1

дА?2

дУ

дУ

ср

д^р

х е [0; /1] , у е [0; /2] , *> 0, ^р(х, 0, *) = ^р(х, /2, *) = ^р(0, у, *) = ^р(/1, у, *) = 0, 1>р(х, у, 0) = 0.

(29)

(30)

д_ дж

+— дж

д_ ду

Е л

дЛт(ТЛ дТ ~ (Т^д^р

дж + Л11 (ТЧ "аж

^ ЗД | + Т12 (т )( £)

ду

дЛт(Т Л дТ + г / д^р Л12 I "^Т + Л12 (Т )

дЛ^2

дж

дж

+

+

+

А

ду

Лр (Т) + £ *^ ду + Т22 (Т>(£)

т=1

= ср

ди>р

д*

ж е [0; /1] , у е [0; /2] , ^р(ж, 0,*) = ^р(ж,/2,*) = и^(0, у, *) = ^р(/1,у,*) = 0, ^р(ж,у, 0) = 0.

Задачи (26)

(32)

(33)

(34) -(4)

8), (29)-(31) и (32)-(34) решаются совместно с прямой задачей ( с помощью метода переменных направлений с экстраполяцией, причём на каждом временном слое используются значения температур, найденные при решении прямой задачи (1)-(4). Таким образом, одновременно находятся элементы матрицы (18) и компоненты вектора (20), которые подставляются в выражение (22).

Фундаментальной проблемой при решении обратных задач является корректность метода решения [1, 5], поскольку коэффициентные обратные задачи зачастую могут быть некорректными, особенно при использовании экспериментальных значений с погрешностями. Здесь можно применить метод регуляризации А.Н. Тихонова [10], однако в этом случае сложно подобрать параметры регуляризации. С другой стороны, если экспериментальные значения считать точными, то можно доказать теоремы о достаточных условиях существования и единственности обратной коэффициентной задачи теплопроводности [4]. Более конструктивным подходом является численный эксперимент, когда при решении обратной задачи к экспериментальным значениям добавляются погрешности до тех величин, пока решение не станет неустойчивым. Именно такой подход и предлагается в настоящей работе.

3. Результаты численных решений

В качестве иллюстрации работоспособности изложенного метода и программного комплекса восстанавливаются нелинейные компоненты тензора теплопроводности углерод-углеродного композитного материала с армированным однонаправленно непрерывным волокном под углом р = 30° к одной из границ пластины. Главные компоненты тензора теплопроводности имели следующие функциональные зависимости:

Л? (Т) = 0.46875 ■ 10-5Т2 — 0.625 ■ 10-2Т + 5.0625,

Лп (Т) = —0.625 ■ 10-5Т2 + 0.0125Т — 4.25.

Итерационный процесс без погрешностей и с добавленной погрешностью в Т^

Номер итерации 5 А11 ^2 ^2 А?1 ^2 А?2 ап

Без погрешности

0 406181 2 2 0 2 2 0 1

1 12342.1 2.84631 1.54525 0.718053 2.98039 2.41681 0.565249 1.2

2 203.876 2.18072 1.55032 0.824741 3.68822 2.44946 0.939521 1.44

3 7.1019 2.42013 1.55114 0.798976 3.44452 2.45603 0.840014 1.728

4 5.96245 2.44688 1.55192 0.802904 3.44453 2.45661 0.850996 2.073

5 5.96186 2.4473 1.55199 0.802965 3.44458 2.45666 0.851008 2.488

6 5.96186 2.4473 1.55199 0.802966 3.44458 2.45666 0.851009 2.985

С добавленной погрешностью

0 412407 2 2 0 2 2 0 1

1 13107 2.82382 1.38153 0.720785 2.85262 2.38604 0.444623 1.2

2 1144.68 1.98963 1.3996 0.820823 3.70417 2.41087 0.944771 1.44

3 595.848 2.28734 1.39834 0.755414 3.40942 2.41461 0.738225 1.728

4 569.87 2.37226 1.39837 0.781108 3.3967 2.40425 0.816026 2.073

5 569.418 2.36953 1.39864 0.779428 3.38991 2.40567 0.809205 2.488

6 569.489 2.37102 1.39862 0.779979 3.39048 2.40538 0.810346 2.985

7 569.489 2.37066 1.39861 0.779769 3.39037 2.40542 0.810189 3.583

Используя соотношения (1), получим функциональные зависимости для компонентов тензора

ЛЦ (Т) = 0.1953125 ■ 10-5Т2 - 0.15625 ■ 10-2Т + 2.734375,

Л22 (Т) = -0.3515625 ■ 10-5Т2 + .78125 ■ 10-2Т - 1.921875,

Л12 (Т) = 0,4736076428 ■ 10-5Т2 - 0,8118988162 ■ 10-2Т + 4.032430788.

Для решения прямой задачи с граничными условиями Т|г = Ттах = 1400 К, Т (ж, у, 0) = Ттт = 600 К принимались входные данные ср = 2.25 ■ 106 Дж/(м3- К), /1 = 0.1 м, ¿2 = 0.06 м. Экспериментальные значения температур вычислялись по этим данным при решении задачи (1)-(4) в точках ж = {0.01; 0.05; 0.09} , у = {0.004; 0.012; 0.02} в моменты времени ^ = 30, 40, 50, 60, 70, 80, 90, 100 с.

Температурный интервал разбивался на три элемента (М = 4): Тт = {600; 900; 1100; 1400}, т = 1,4. Таким образом, определялись 12 искомых параметров: Л™, Лт>, Лт>, т = 1, 4. В таблице приведены итерационные процессы по восстановлению коэффициентов Л11, Л*2, Л12, Л31, Л22, Л?2 без погрешностей в экспериментальных данных Т^

АТ,-

< 5 с равномерным

и с добавленной абсолютной погрешностью в интервале распределением. Видно, что в обоих случаях функционал стремится к стационарному значению. Начальные величины искомых параметров значительно (в 2-3 раза) отличаются от искомых, и тем не менее итерационный процесс сходится с приемлемой скоростью (за шесть — семь итераций), причём параметры на нулевой итерации выбраны из предположения, что материал исследуемой пластины является изотропным с теплопроводностью Л = 2 Вт/(м ■ К).

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

600 900 1100 1400 600 900 1100 1400

в

Рис. 2. Восстановленные компоненты Л11 (Т) (а), Л22 (Т) (б), Л12 (Т) (в) тензора теплопроводности при наличии погрешностей в определении экспериментальных значений Т *

с учётом погрешностей в экспериментальных значениях температур Т * (кривая 3 — АТ,* е [—5К; 5К], кривая 4 — АТ,* е [—10К; 10 К]) и нелинейных коэффициентов, заложенных в экспериментальные значения ТД- (кривая 1). Из рисунков видно, что несмотря на добавленную погрешность в экспериментальные данные восстановленные компоненты тензора теплопроводности достаточно точно аппроксимируют кривые, соответствующие нелинейным коэффициентам, заложенным в экспериментальные значения Т*.

Выводы

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

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

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

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

[1] Алифанов О.М. Идентификация процессов теплообмена летательных аппаратов (введение в теорию обратных задач теплообмена). М.: Машиностроение, 1979. 216 с.

[2] Бек Дж., Блакуэлл Б., Сент-Клэр Ч., мл. Некорректные обратные задачи теплопроводности: Пер. с англ. М.: Мир, 1989. 312 с.

[3] Самарский А.А., Вабищевич П.Н. Численные методы решения обратных задач математической физики. М.: Изд-во ЛКИ, 2009. 480 с.

[4] Кузнецова Е.Л. Восстановление характеристик тензора теплопроводности на основе аналитического решения задачи теплопереноса в анизотропном полупространстве // Теплофизика высоких температур. 2011. Т. 49, № 6. С. 1-8.

[5] Формалёв В.Ф. Тепломассоперенос в анизотропных телах. Обзор // Там же. 2001. Т. 39, № 5. С. 810-832.

[6] Криксин Ю., Плющев С.Н., Самарская Е.А., Тишкин В.Ф. Обратная задача восстановления источника для уравнения конвективной диффузии // Математическое моделирование. 1995. Т. 7, № 11. С. 95-108.

[7] Формалёв В.Ф., РЕвизников Д.Л. Численные методы. М.: Физматлит, 2004. 400 с.

[8] Формалёв В.Ф., Колесник С.А., ЧипАшвилли А.А. Численное моделирование теплопереноса в анизотропных телах с разрывными характеристиками // Там же. 2004. Т. 16, № 5. С. 94-102.

[9] Тихонов А.Н., Леонов А.С., ЯголА А.Г. Нелинейные некорректные задачи. М.: Наука, 1995.

Поступила в 'редакцию 14 августа 2012 г., с доработки — 28 декабря 2012 г.