Метод Галеркина для решения задачи теплопереноса при исследовании условий возникновения критических тепловых режимов при ламинарном движении химически реагирующей среды

ПРОБЛЕМЫ ЕСТЕСТВЕННЫХ НАУК

УДК 532.135

А.И. Кадыйров, Е.К. Вачагина

МЕТОД ГАЛЕРКИНА ДЛЯ РЕШЕНИЯ ЗАДАЧИ ТЕПЛОПЕРЕНОСА ПРИ ИССЛЕДОВАНИИ УСЛОВИЙ ВОЗНИКНОВЕНИЯ КРИТИЧЕСКИХ ТЕПЛОВЫХ РЕЖИМОВ ПРИ ЛАМИНАРНОМ ДВИЖЕНИИ ХИМИЧЕСКИ РЕАГИРУЮЩЕЙ СРЕДЫ

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

Метод Галеркина, ламинарное течение, химическая реакция, критический режим.

АХ Kadiyrov, E.K. Vachagina

GALERKIN METHOD FOR SOLVING THE PROBLEM OF HEAT TRANSFER CONDITIONS IN THE STUDY OF CRITICAL THERMAL RATES IN THE LAMINAR FLOW OF CHEMICALLY REACTING LIQUID

This paper presents Galerkin method for solving the problem of finding critical conditions of thermal rates in the laminar flow of chemically reacting non-Newtonian liquid. Numerical calculations show that the transition goes in a very narrow interval of changing the parameter characterizing the initial state of the system from slow mode to explosive. The boundary intervals of safe thermal modes of heat transfer are presented in the article.

Galerkin method, laminar flow, chemical reaction, a critical mode.

Введение

Прикладное значение математической теории «теплового взрыва» чрезвычайно велико для безопасной эксплуатации теплотехнологического оборудования, использующего в качестве рабочих сред химически реагирующие среды. Основы этой теории были заложены в трудах Н.Н. Семенова, Д.А. Франк-Каменецкого, Я.Б. Зельдовича, Г.И. Баренблатта,

О.М. Тодеса [1-3].

Исследования процессов теплового взрыва имеют важное значение для решения как научных, так и практических задач. Эти исследования получили широкое развитие и отражены в монографиях и обзорных статьях, например, [4]. Одной из задач о тепловом

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

По настоящее время явление теплового взрыва вызывает интерес в научной сфере, подтверждением служат публикуемые в научных журналах статьи по данной тематике [510].

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

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

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

% %+* %=2 А (*, е)% )-|+| (*, е)%)+

% % +т % =1!¿(е)Ш+4¿(е)Ш+2д(¿(е)Ш -

х дх г дг дх ^ 2 ’ дг ) дх ^ 2 ' дх ) дг ^ 2 ’ дг )

—+1 ¿

дг х У 2 ! дх х К ! дх

А™-1+^ !)=Щ|>8 ех1)1+зё+^

= 0.

х дх дг

Краевые условия:

Условие прилипания жидкости при х = 1:

Жг = 0, Жх = 0.

Условие на оси трубы при х = 0:

Щ = 0, = 0, дде = °-

дх дх

Условие на входе в трубу при г = 0 :

е = е0, Щ Щ 0

где е0 - безразмерная температура на входе, 0

с Щ = 0.

Условие для температуры при х = 1:

д 0 *

^ = -Б1 •

дх

Здесь х, г - координаты; Жх, - радиальная и осевая составляющие вектора

скорости; д - безразмерная динамическая вязкость жидкости; 12 - безразмерный второй

дР

инвариант тензора скоростей деформации;------------безразмерный перепад давления; 5 -

дх

безразмерный параметр, характеризующий тепловыделение от химической реакции; х -параметр, характеризующий диссипативный источник тепловыделения; а - параметр, характеризующий отношение энергии активации вязкого течения к энергии активации химической реакции; в - параметр, характеризующий чувствительность системы к окружающей среде; у - параметр, характеризующий отношение расхода жидкости к длине дуги канала, помноженный на коэффициент температуропроводности; Е, - безразмерный параметр; Ы - число Био; с0, с1 - реологические параметры, характеризующие свойство *

жидкости, е - температура окружающей среды.

Метод решения

Рассмотрим метод решения системы уравнений (1)-(4) с краевыми условиями (5)-(8) для высоковязких сред, для которых числа Рг >> 1. Для выбранных сред время гидродинамической релаксации много меньше тепловой, поэтому процесс можно считать квазистационарным, соответственно профиль скорости почти мгновенно подстраивается под профиль температуры. Кроме того, изменение соответствующих составляющих вектора скорости в поперечных направлениях много больше их изменений в направлении основного движения.

В силу того, что рассматриваемая система уравнений (1)-(4) с краевыми условиями (5)-(8) носит нелинейный характер, то аналитически найти решение не представляется возможным. В качестве численного метода решения предлагается рассмотреть применение метода Галеркина. Численная реализация решения поставленной задачи будет носить итерационный характер, смысл которого заключается в разбиении оси г на слои с заданным шагом.

о, = {гк|гк = ¿‘-'+к‘, к = Щ} ,

где г0 = 0; Ик = 0,01 - шаг сетки разбиения оси г ; - выбирается в зависимости от

длины рассматриваемого участка.

Таким образом, на каждом слое сетки поочередно решались уравнение переноса энергии (3) и упрощенное уравнение движения вида

1 дх[х¿

х дх ^ дх ) дг 0 )

Уравнение неразрывности (4) используется для определения радиальной

дР

составляющей Жх, а уравнение (1) - для определения — на каждом слое по г .

дх

В качестве реологического уравнения состояния среды принято уравнение Кутателадзе - Хабахпашевой [11], которое в безразмерном виде имеет вид

д

1+ P0J С - (с0 -1) • exp (- с • д(/2 , 0))7Т) ’

где с 0, с 1 - реологические параметры, характеризующие свойства жидкостей.

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

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

Далее производится расчет гидродинамических характеристик потока на этом же слое на основании найденного поля температур 0, принимая в первом приближении значение dP/ dz равным значению dP/ dz на предыдущем слое. После расчета гидродинамических полей определяется уточненное значение dP / dz и распределение вязкости.

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

Для решения уравнения (10) применяется итерационный метод решения. При этом на нулевом шаге итерации (j = 0) считаются заданными значение температуры, а профиль скорости является решением стационарной задачи (10), (5), (6) с Wx = 0. Смысл итерационного метода заключается в следующем. Вначале определяется интервал изменения A <dP / dz < B . Затем применяем метод, аналогичный итерационному методу «деления отрезка пополам», а именно, уравнение (10) решается с dP / dz = (A + B ) /2. Полученное решение WzJ используем для пересчета dP / dz. Решение повторяется и т.д. Процесс продолжается до тех пор, пока:

WJ - W/

WJ

< s процесс сошелся, получено искомое

решение.

Если процесс завершается на этапах 1) или 2), то необходимо изменить значения концов интервала (A, B), применяя метод, аналогичный итерационному методу «деления отрезка пополам». Пересчитываем значения концов интервала C = (A + B)/2 . Если при

dP/ dz = C выполняется первое условие 1), то A = C; если условие 2), то B = C; и т.д.

Решение уравнения движения. Рассмотрим подробнее решение уравнения (10) с dP/dz = const. Вследствие того, что уравнение (10) является нелинейным, для его решения используется численный итерационный алгоритм, который заключается в применении метода Галеркина к уравнению (10) и решении получившейся системы алгебраических уравнений методом Гаусса.

Уравнение (10) в операторной форме представимо как

L

где

д ( ди} ци) = дХ ^ к (х)^ии ) + Р (х,и).

В частном случае имеем

и = ^, К(х) = хд, Р(х,и) = х • |^Гии'1 хёх.

0 ^дх )

Для решения уравнения типа

X = хг-1 + К,г = 1, , К = N!

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

Решение ищется методом конечных элементов, где в качестве базисных функций используются кубические сплайны [12]

Фг (X) = ф ^ I, М(х) = КМ ^ I , г = 1, ... , N -1.

К

ф0(х) = тах1 0, Кф

V0(х) = тах1 0, К у

х - а

~К~ ,

' х - а'

~К~.

К

, МN (х) = тах І0, К ф

VN (х) = тах1 0, К у

Ь - х '

. К , Ґ Ь - х'

~К~.

ф(і) =

- 2(і +1)3 + 3( +1)2, і є [-1, 0],

2(і -1)3 + 3(і -1)2, і є [0, 1], у(і) =

0,

і £ [-1, 1]

(і +1) - (і +1) , і є [-1, 0], (і -1)3 + (і -1)2, і є [0, 1]

0,

і £ [-1, 1]

Как следствие, всякая функция ик (х) е У2 м ^ ^[а, Ь] может быть единственным образом записана в виде

N

и * иЬ (х) = 2 (сгфг (х) + С>г (х)),

г=0

где коэффициенты сг суть значения функции ик (х) в узлах сетки (15), а с' - значения ее

первой производной на сетке (15).

Подставив вместо искомой функции ее разложение вида (20), получим дискретный оператор

Ґ

К д

ь(ип) ^

дх

N

д1 !

К ( х)А^-----------------------

дх

Л

+ Е (х,иК).

Воспользуемся условием ортогональности невязки базисным функциям (Ь(ик), ф; (х)) = 0, (Ь(ик),у;(х))= 0 ] = 0,1,...,N.

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

(ь(ик X V}.(х) )= | Ь(иН) • V}.(х) ёх, і = ^..^ N.

Таким образом, учитывая (21), (19), решение операторного уравнения вида (12), (13) состоит в решении следующей системы уравнений

СО = 1х

Е'c fj діK (x) дфд(x)1 •ф j(x) dx Vic;f k (x) dyd(x) 1 •ф j(x) dx

f 0

dx I

dx

f 0

dx I

dx

= -J F (x,Uh) •ф j (x) dx, j = 1,..., N -1;

Ec fJ dx fK (x) дф^ j(x) dx l+Ec;fJ dx fK (x) ^ )^ j(x) ^

¿=0 f 0 dx f dx ) ) i=0 f 0 dx f dx

)

f0

= -J F ( x, Uh) • у j ( x) dx .

)

Рассмотрим численную реализацию решения уравнения (10) с краевыми условиями (5), (6) на фиксированном слое ¿к, к > 0, к є N. После упрощения системы уравнений (24), (25) с учетом (13) получим следующую систему алгебраических уравнений относительно

неизвестных Cj , С j

aj-1cj-1 + ajcj + aj+1cj+1 + bj -1cj-1 + bjC j + bj+1cj+1 = -f

q.-Cj-1 + . + qj+1'j.+1 + ij.-'-1 + Sj'j + sj+1cj+1 = - gj

(j = 1, 2,..., #-1).

c0 c h , . , n

- — + — +—c0 + 4c1 = 0,

2 2 10 0 1

11h2 h h2 , n

------c0 +------c----c = 0 .

210 0 10 1 60 1

9h 70 13h2 420

c,r, +

13h

35

’cN +

h2 420

cN-1 -'

11h2

210

c'n = 0,

c ---------

N-1

11h2 210

' cN -

h3

140

cN-1 +

h3

105

cN = 0.

Получившуюся систему алгебраических уравнений решаем методом Гаусса. Для вычисления интегральных выражений коэффициентов системы (26)-(27) используется квадратурная формула Симпсона.

Решение уравнения неразрывности. Рассмотрим теперь решение уравнения неразрывности (4) на фиксированном слое Zk, к > 0, k е N. В качестве метода решения выбран метод Галеркина. Сначала непрерывные функции аппроксимируем дискретными на сетке (14): W, ~(W,), , W, ~ (W,),, при этом x ~ xi. При численной реализации решения

данного уравнения значение второго слагаемого (W,) берется из предыдущих двух слоев,

dW- dW.

поэтому на фиксированном слое ------------ = const. Отсюда можно записать --------- = M ,

dZ dZ

M| < +да . С учетом этого уравнение (4) возможно записать в операторном виде (12) для

U = Wx.

d

L(U) = —(xU) + xM .

dx

Применив метод Галеркина и учитывая краевые условия (5), (6), получим следующую систему алгебраических уравнений относительно неизвестных с, с'

((-27 + 9^_j)cj.-1 + 26c. - (27 + 91 j).)+

h

i=0

i=0

где М =

+

-3 ((-\0 + 21/) ( + \4/С/ + (\0 + 21/) /)+

Их

+

13И 9И \\И2 , \3И2 , л

-----Со + —С\ +---с0--------------с\ = 0,

35 0 70 1 210 0 420 1

пи2 ш2 и3 , И , = п

С0 +------------с +--С0-----------с = 0,

210 0 420 1 105 0 140 1

9И 13И И2 , 11^2 ,

—с*-1 +---------с* +-------слг-1---------с* = 0,

70 *1 35 * 420 210 *

13И2 11И2 И3 , И3 , л

--------с*-1------------с*-------слг-1 +----с* = 0,

420 210 * 140 105 *

( дЖ

1 = ((Ж/-(Ж, У/Г1), (/ — 1, 2, *-1), (I — 2, 3, Л).

Систему (32)-(38) решаем методом Гаусса.

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

ик (х, г) е У2* с помощью рядов вида

*

и « иИ

1=0

здесь коэффициенты с1 (,) суть значения функции ик (х, г) в узлах сетки (15), а с'(,) -

значения ее первой производной по х на сетке (15). Здесь в качестве базисных функций взяты кубические сплайны (19).

В операторной форме уравнение переноса энергии можно записать в виде (11),

где

ди д ( ди 1

ци) - К0

дг дх V дху

К0

Р

дх 1 + ри

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

_м ( дс ( ,) 1 дсг ( ,) 1 1

Ё —ь-1К0(х)фг(х)ф;(х)дх+—д— 1К0(х)^-(х)ф/(х)дх I-

1=0 V дг 0 дг 0 у

- ё(с(г )| дх {К1(х) дфд( х) 1 ф /(х) дх+1 с'( г) дх(К1(х) дфдхх) И ф /(х) дх=

¿=0 V 0 дх V дх у 0 дх V дх уу

= 1Р

Ё

1=0

| К„ (х) ф, (х) х|/ (х) ¿/х + 1 Кп (

02 п 02 п 6

- Е [с, к1 (х)%~] (х)й6г+с>^ 4

2=о V о дх\ дх ) ,

= | Р(х, иь, ) ф; (х) ¿/х,

0

где 7 = 0,...,#.

С учетом краевых условий (8), (6), получаем следующую систему обыкновенных дифференциальных уравнений

относительно неизвестных с

7

Зсн

-/ + ~ЗГ“/ +

дг

дг

дг

дг

0-1

дг

дг

+^3 - Т />1 +(-Т - I(5+^ +

+10И ( -29 /у / - Ис; +10И (+29/) = - /1,

дс

—------Ян +—/—Я/ +-/-----/ + —-/ + ——^ + —-------

дг

дг / дг / 1 дг /-1 дг

+ ЛИ(-7 + 9/)с/-1 -ТС ^И(7 + 9/)с-+' +

дг

10

10

7+1

+30 И'(-15+16/И--15 и2-/'у; - 30 И2 (5+^-/^ = -^/

(/ = 1,..., * -1),

£°(£) + уу; (г)+4с; (г) = 0,

2 10 1

2

11И2 210

И И2

с0( г) + ТТ с1( г) - с; (г) = 0 =

10

( 1 9Б1*И 1

— + -2

70

с*-1( г

60

(1 13Б1*И 1

- + -

35

с* ( г ) -

И

(- 42 +13 Б1* И) 1 (г) + — (21 - ПБ^иС (г) = - 0*, 42^ ' ^ 210'' ' 2

13

420 ' 21^ ' *

— (1 + Бг*) с* 1(г) + — (1 + Бг*) с* (г) = -—0*. 140 у ’ 105 у ’ 12

Полученная система обыкновенных дифференциальных уравнений (45)-(50) решена с помощью метода Рунге - Кутта с условием (7).

2

Рис. 1. Зависимость безразмерной температуры от параметра Франк-Каменецкого; аналитический метод - 1; метод Галеркина - 2

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

предложенной Франк-Каменецким в работе [1]. На рис. 1 представлены графики аналитического и численного решения задачи о самовоспламенении.

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

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

Результаты численных исследований

Основная задача математической теории теплового взрыва заключается в исследовании динамики исследуемого процесса при заданных значениях управляющих параметров. Для классической модели теплового взрыва эти характеристики отражает некоторый параметр (например, 5 -параметр Франк-Каменецкого), значение которого определяется начальным

состоянием химической системы. В

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

Численные расчеты по исследованию теплопереноса при ламинарном движении химической реагирующей высоковязкой

среды по модели (1)-(8) показывают, что переход от медленного режима к взрывному происходит в чрезвычайно узком

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

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

Границы интервалов безопасных тепловых режимов теплопереноса

Рис. 2. Распределение безразмерной температуры по длине канала на начальном тепловом участке

I - VIII режимы

I II III IV V VI VII VIII

5 < 1,1 0,98 1,08 1 1,1 0,67 1,09 1,07

X < 0,01 0,1 0,01 0,01 0,01 0,01 0,01 0,01

а > 1,27 1,27 0,63 1,27 1,27 1,27 1,27 1,27

e > 0,156 0,156 0,156 0,1 0,156 0,156 0,156 0,156

Y > 0,03 0,03 0,03 0,03 0,1 0,03 0,03 0,03

Bi* > 2 2 2 2 2 1 2 2

Co 3 3 3 3 3 3 1 0,4

Ci 0,04 0,01 0,04 0,04 0,04 0,04 1 3

Выводы

1. Представлен алгоритм решения задачи теплопереноса при исследовании условий возникновения критических тепловых режимов при ламинарном движении химически реагирующей среды на основе метода Галеркина.

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

Работа выполнена при проведении госконтракта № 1115 от 26.08.2009 в рамках реализации ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009-2013 годы, а также при финансовой поддержке РФФИ (№ 08-08-12109 офи).

ЛИТЕРАТУРА

1. Франк-Каменецкий Д.А. Диффузия и теплопередача в химической кинетике / Д. А. Франк-Каменецкий. М.: АН СССР, 1947. 492 с.

2. Конторова Т.А. К теории теплового взрыва / Т.А. Конторова, О.М. Тодес // Журнал физической химии. 1933. Т. 4. С. 81.

3. Тодес О.М. Теория теплового взрыва / О.М. Тодес // Журнал физической химии. 1939. Т. 13. С. 75.

4. Мержанов А.Г. Современное состояние теории теплового взрыва / А.Г. Мержанов, Ф.И. Дубовицкий // Успехи химии. 1966. Т. 36. Вып. 4. С. 656-683.

5. Новожилов Б.В. Тепловой взрыв при ламинарном движении реагирующей среды / Б.В. Новожилов, Н.Г. Самойленко, Г.Б. Манелис // Химическая физика. 2004. Т. 23. № 3. С. 62-66.

6. Штейнберг А.С. Тепловой взрыв при наличии колебаний температуры среды, окружающей реакционную смесь / А.С. Штейнберг, С.И. Худяев // Доклады РАН. 2005. Т. 403. № 2. С. 212-215.

7. Новожилов Б.В. Тепловой взрыв в перемешиваемой среде / Б.В. Новожилов, Н.Г. Самойленко, Г.Б. Манелис // Доклады РАН. 2002. Т. 385. № 2. С. 217-219.

8. Шкадинский К.Г. Анализ теплового взрыва в системах «пористый реагент -активный газ - твердый продукт» / К.Г. Шкадинский, Н.И. Озерковская, С.В. Костин // Физика горения и взрыва. 2004. Т. 40. № 3. С. 40-50.

9. Кравченко В.Ф. Асимптотики Пуанкаре решений задач нерегулярного тепло- и массопереноса / В.Ф. Кравченко, Г. А. Несененко, В.И. Пустовойт. М.: Физматлит, 2006. 216 с.

10. Котович А.В. Решение многомерной нелинейной задачи об очаговом тепловом взрыве асимптотическим методом / А.В. Котович, Г.А. Несененко // Материалы VI Минского Междунар. форума по тепломассообмену. Минск, 2008. С. 268-269.

11. Кутателадзе С.С. К гидродинамике жидкостей с переменной вязкостью / С.С. Кутателадзе, В.И. Попов, Е.М. Хабахпашева // Прикладная математика и техническая физика. 1966. № 1. С. 45-49.

12. Бахвалов Н.С. Численные методы / Н.С. Бахвалов. М.: Наука, 1975. 632 с.

Кадыйров Айдар Ильдусович - Kadiyrov Aydar Ildusovich -

кандидат технических наук, научный Candidate of Technical Sciences,

сотрудник Исследовательского центра Member of the Research Center

проблем энергетики Учреждения РАН of Power Engineering Problems

Казанского научного центра РАН

Вачагина Екатерина Константиновна -

доктор технических наук, заведующая лабораторией теплофизических исследований

Исследовательского центра проблем энергетики Учреждения РАН Казанского научного центра РАН

of Russian Academy Institution of Kazan Scientific Centre RSA

Vachagina Ekaterina Konstantinovna -

Doctor of Technical Sciences,

Head of Heat-transfer research Laboratory of the Research Center of Power Engineering Problems of Russian Sciences Academy of Kazan Scientific Centre RSA

Статья поступила в редакцию 03.12.09, принята к опубликованию 30.06.10