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

МЕТОДЫ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ

УДК 536. 2 (075)

МЕТОД ДОПОЛНИТЕЛЬНЫХ ИСКОМЫХ ФУНКЦИЙ В ЗАДАЧАХ ТЕПЛОПРОВОДНОСТИ С ПЕРЕМЕННЫМИ ФИЗИЧЕСКИМИ

СВОЙСТВАМИ СРЕДЫ1

Е.В. КОТОВА, А.В. ЕРЕМИН, В.А. КУДИНОВ, В.К. ТКАЧЕВ, А.Э. КУЗНЕЦОВА ФГБОУВО «Самарский государственный технический университет», г. Самара, Российская Федерация E-mail: [email protected]

Авторское резюме

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

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

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

1 Работа выполнена при поддержке РФФИ в рамках научного проекта №18-38-00029 мол._а.

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

METHOD OF ADDITIONAL UNKNOWN FUNCTIONS IN HEAT CONDUCTIVITY PROBLEMS WITH VARIABLE PHYSICAL PROPERTIES OF THE MEDIUM

E.V. KOTOVA, A.V. EREMIN, V.A. KUDINOV, V.K. TKACHEV, A.E. KUZNETSOVA Samara State Technical University, Samara, Russian Federation E-mail: [email protected]

Abstract

Background. Finding analytical solutions to the problems of thermal conductivity with variable physical properties of the medium by classical analytical methods is very complicated mathematically. The known expressions representing complex infinite series including two types of Bessel functions and gamma-functions are, in fact, numerical as they require a numerical solution to complex transcendental equations with eigenvalues of the boundary problem. Such solutions can hardly be used in engineering applications, especially in cases when a solution to a certain problem is only an intermediate stage in other problems (such as thermoelasticity and control problems, inverse problems, etc.) which can be solved effectively only by finding analytical solutions to the initial problems. Therefore, an urgent problem now is to develop new methods of obtaining analytical solutions to the abovementioned problems, at least approximate ones. Materials and methods. The study employed methods of additional boundary conditions and additional unknown functions in the integral method of heat balance.

Results. High-precision approximate analytical solutions to the transient heat conduction problem with non-homogeneous physical properties of the medium for an infinite plate under symmetric boundary conditions of the first type have been obtained. The initial problem for partial differential equations is reduced to two problems in which ordinary differential equations are integrated. Additional boundary conditions are defined in such a way that their fulfillment in accordance with the new method is equivalent to the result of solving the initial partial differential equation at the boundary points and at the temperature perturbation front (for the first stage of the process).

Conclusions. By combining methods with finite and infinite heat propagation rate we have been able to obtain high-precision analytical solutions for the whole time range of the unsteady process including its small and ultra small values. The solutions look like simple algebraical polynomials not including special functions (Bessel, Legendre, gamma-functions and others). Since it is not necessary to directly integrate the initial equations by the space variable and to reduce them to ordinary differential equations with additional unknown functions, the considered method can be used for solving complex boundary problems in which differential equations do not allow distinguishing between the variables (into nonlinear, with linear boundary conditions and heat sources, etc.).

Key words: transient heat conduction, variable physical properties of the medium, additional unknown functions, additional boundary conditions, temperature perturbation front

DOI: 10.17588/2072-2672.2019.2.059-070

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

Известные решения такого типа задач, найденные классическими методами,

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

емые аналитические решения являются по сути численными.

Методы исследования. Для нахождения приближенных и точных аналитических решений краевых задач большое распространение получили ортогональные методы взвешенных невязок (Л.В. Канторовича, Бубнова - Галеркина, интегральный метод теплового баланса и др.) [2-10]. Среди них наиболее простым и универсальным является интегральный метод. При его использовании процесс теплопроводности разделяется на две стадии, первая из которых характеризуется постепенным продвижением фронта температурного возмущения от поверхности к центру, а вторая - изменением температуры по всему объему тела. Для каждой стадии вводятся дополнительные искомые функции. Дополнительная функция для первой стадии описывает закономерность перемещения фронта возмущения по координате во времени, а во второй - закон изменения температуры во времени в центре пластины. Таким образом, обе эти функции зависят лишь от времени, и поэтому их введение позволяет свести решение исходного уравнения в частных производных к интегрированию двух обыкновенных дифференциальных уравнений относительно дополнительных искомых функций. Важное преимущество метода в том, что при его использовании практически не накладывается каких-либо ограничений на вид дифференциальных операторов краевых задач - они могут быть нелинейными, с переменными физическими свойствами среды и др. Используя эти методы для решения многих сложных задач математической физики, не поддающихся решению классическими аналитическими методами, в ряде случаев удается получать аналитические решения [5-10].

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

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

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

с( х )у( х)

дТ (x,t) dt

д_ дх

Ц( х)

дТ (x,t) дх

(1)

(t > 0; 0 < х <5);

Т (х, 0) = То; (2)

дТ (0,t)/ дх = 0; (3)

Т (5,t) = Ть (4)

где Т - температура; х - координата; t -время; Ц(х) - коэффициент теплопроводности; с(х) - теплоемкость; у(х) - плотность; Т0 - начальная температура; Т -температура стенки при х = 8; 8 - половина толщины пластины.

Найдем решение задачи (1)-(4) в случае, когда произведение су = const, а коэффициент теплопроводности Ц является экспоненциальной функцией координаты х:

Ц( х) = Ц0 exp(-mx), (5)

где m > 0 - коэффициент, характеризующий интенсивность изменения коэффициента теплопроводности по координате х; Ц0 = const - коэффициент теплопроводности пластины при х = 0.

Введем следующие безразмерные переменные и параметры:

© =Т-Т° ; Fo = a4 ; а0=Ц° ; х ■ (6)

Т - Т0 82 с у 8

где © - безразмерная температура; Fc -число Фурье; а0 - коэффициент темпера-

туропроводности пластины при х = 0; £ -безразмерная координата.

С учетом обозначений (6) задача (1)-(4) принимает следующий вид:

s©(ç,Fo) _ д

д©(£, Fo)

(7)

(8) (9) (10)

5Ро

(Ро > 0; 0 <£< 1) ©(£,0) = 0; 5©(0,Ро)/а£ = 0; ©(1,Ро) = 1 , где V = т8 .

Разделим процесс теплообмена на две стадии по времени: 0 < Ро < Ро и РО < Ро < да . Для этого введем движущуюся во времени границу (фронт температурного возмущения), разделяющую исходную область 0 < £ < 1 на две подобласти щ(Ро) <£< 1 и 0 <£<щ(Ро), где щ(Ро) -функция, определяющая продвижение границы раздела во времени (рис. 1,а).

1,0

©

4 1,0

а)

1,0 ©

1,0

б)

Рис. 1. Расчетная схема теплообмена

Первая стадия заканчивается при достижении движущейся границей центра пластины (£ = 0), т. е. при Ро= Ро. При этом понятие фронта возмущения теряет смысл и начинается вторая стадия теплообмена, в которой изменение температуры

происходит по всей толщине пластины: 0 < £ < 1. Математическая постановка задачи для этой стадии будет рассмотрена ниже. Применительно к первой стадии процесса математическая постановка задачи после введения независимой переменной р = 1 - £ записывается в следующем виде (рис. 1,б):

д©(р, Fo)

dFo

д_

др

(11)

(12)

(13)

(14)

v(i-p) д©(р, Fo) др _

(0 <Fo< Foi ; 0<р<qi(Fo)); © (0, Fo) = 1 ; © (9i, Fo) = 0 ; д©(^ ,Fo)/ др = 0.

Задача (11)-(14) не одержит начальнoгo услoвия (8), так как за пределами фрoнта температурнoгo вoзмущения oна вooбще не oпределена. Так как в начальный мoмент времени ( Fo = 0 ) величина фрoнта температурнoгo вoзмущения равна нулю, то в качестве начальнoгo услoвия принимается услoвие вида q (0) = 0. Задача (11 )—(14) не одержит также граничнoгo услoвия вида (9), так как oнo не влияет на прoцесс теплooбмена в первoй егo стадии.

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

Решение задачи (11)-(14) разыскивается в виде

©(р,Fo) = £ ak (qi) рк ,

(15)

k=0

где ак(щ) - неизвестные коэффициенты, определяемые из условий (12)-(14). После подстановки (15) (при п = 2) в (12)-(14) относительно ак (щ), (к = 0,1,2) получаем систему трех алгебраических линейных уравнений, из решения которой находим а0 = 1; а1 = -2/д, ; а2 = 1/щ2 .

С учетом найденных значений коэффициентов ак (щ), (к = 0,1,2) соотношение (15) примет вид

©(р^о) = (1-р/Я1? . (16)

Для нахождения неизвестной функции щ(Ро) составляется невязка уравне-

0

0

р

ния (11) и определяется интеграл от нее в пределах толщины прогретого слоя (интеграл теплового баланса):

91

д0(р^о) -Ро

91 я

Ф=[ -

J Лр

у(1-р) д©(р,Ро) ф

dр. (17)

0 - " 0

Вычисляя интегралы в (17), получаем

= 6exp(-v)dFо . (18)

Решение уравнения (18) при начальном условии д (0) = 0 имеет вид

д,(Ро) = 2Л/3роёхру /ехрV . (19)

Положив д (Ро^ = 1 (при V = 0,01), находим время Ро = 0,084, при котором фронт температурного возмущения достигает координаты р = 1.

Формулы (16), (19) являются решением задачи (11)-(14) в первом приближении. Сравнение результатов расчетов по формуле (16) с решением методом конечных разностей позволяет заключить, что в диапазоне 0,025 < Ро < 0,005 их расхождение не превышает 6 %. Повышение точности связано с увеличением числа членов ряда (15), неизвестные коэффициенты которого будем находить из основных (12)-(14) и дополнительных граничных условий. Отметим, что во втором и каждом последующем приближении необходимо использовать по три дополнительных граничных условия - одно при р = 0 и два других при р = д(Ро). Для их получения будем последовательно дифференцировать граничные условия (12)-(14) по переменной Ро, а уравнение (11) по переменной р. Сравнивая получающиеся при этом соотношения, можно найти необходимое количество дополнительных граничных условий. Для получения первого из них продифференцируем (12) по Ро: 5©(0^о)

ЛРо

= 0.

(20)

Уравнение (11) запишем для точки

р = 0:

50(0,Ро) ЛРо

= ve

5©(р, Ро)

Ф

+ е

52©(р,Ро)

р=0 ф р=0 (21)

Сравнивая (20) и (21), получаем первое дополнительное граничное условие:

д2©(0,Ро) 50(0,Ро)

ф2

Ф

Для получения второго дополнительного граничного условия продифференцируем соотношение (13) по переменной Ро:

5©(р,Ро)

5Ро

= 0 .

р=91

Запишем

(23)

уравнение (11) для

р = 9 (Ро): 5©(р,Ро)

5Ро

= ve

р=91

v(1-g1) д©(р,РР) ф

р=91

+е

v(1-q1) 52©(р,Ро) ф2

(24)

р=91

Сравнивая (23) и (24), с учетом (14) получим второе дополнительное граничное условие:

52©(91,Ро)

Ф2

= 0 .

(25)

Для получения третьего дополнительного граничного условия продифференцируем соотношение (14) по переменной Ро:

52©(р,Ро)

ЗРоф

= 0 .

(26)

р=91

Дифференцируя уравнение (11) по переменной р и применяя полученное соотношение к точке р=д(Ро), с учетом (14) и (25) получаем

Л2 ©(р,Ро)

дРоф

р=91

©(р,Ро)

" е я 3

Ф

(27)

р=91

Сравнивая (26) и (27), находим третье дополнительное граничное условие:

53©(91,Ро)

Ф3

= 0 .

(28)

Подставляя (15), ограничиваясь шестью членами ряда, в основные (12)-(14) и дополнительные (22),(25),(28) граничные условия относительно ак(д), (к = 0,5) будем иметь систему шести алгебраических линейных уравнений:

(а0 + а,р + а2р2 + азр3 + а4р4 + а5р5) = 1,

\ /р=0

а0 + а1д^1 + а2^2 + а3д13 + а4д14 + а5д15 = 0, а1 + 2а2^ + 3а3^!2 + 4а4д!3 + 5а5д14 = 0, I (29) а1v- 2а2 = 0,

а2 + 3а3д1 + 6а4^2 + 10а5д13 = 0,

а391 + 6а4Я1 + 10а5^2 = 0.

Ввиду высоких степеней производных в дополнительных граничных услови-

+

ях система уравнений (29) имеет цепочный вид, что существенно упрощает процесс получения ее решения. Из решения системы уравнений (29) находим:

а0 = 1; а = 20 /(дг); а2 = а /2; аз = 20 (уд + 2)/( д3г); а4 = 5(3уд + 8)/(д4 г); а5 = 4(уд + 3)/д5г; г = уд + 8.

Подставляя (30) в (15), получаем

(30)

®(р,Ро) =

(12р + 8д + уд + 4уд1р)(р - д) ^ (уд + 8)д5 '

Подставляя (31) в интеграл теплового баланса (17), относительно неизвестной функции д(Ро) получаем следующее нелинейное обыкновенное дифференциальное уравнение: с^ = 60уд + 480 dРо у2д3 +1 буд2 + 48д '

Следуя методу [13], находится приближенное аналитическое решение уравнения (32):

д#о) = , (33)

А "

2

где ц = 1 (Ро* )* ; Р01

у

90

1280

1. 40 ;

Х = А -

209у3 + 1312у2 + 2304 192у2 + 3072у + 921б

А =

Г81у5 + 4104у4 + 83008у3 + 80537бу2 3317б00у + 265420800

3б21888у + 597196

12

3317600у+2б5420800у

Например, при у = 0,01 ц = 4,47235; Х = 0,5; Ро^ = 0,05 - время достижения фронтом температурного возмущения координаты р = 1 во втором приближении.

Результаты расчетов перемещения фронта температурного возмущения по формуле (19) (решение уравнения (18)), а также по формуле (33) (решение уравнения (32)) в сравнении с решением уравнения (32) численным методом, приведены на рис. 2. Их анализ показывает практическое совпадение численного решения с результатами, полученными по формуле (33). Следовательно, она с достаточной степенью точности может быть использована в качестве аналитического решения уравнения (32).

1,0

р

0,8 0,6 0,4 0,2

2

/ о 1

/ /и /

// о/

ро*=0,05 Ро1=0,084-^

0 0,014 0,028 0,042 0,056 Ро 0,084

Рис. 2. Перемещение фронта температурного возмущения д(Ро) по координате р во времени Ро (у = 0,01): 1 - по формуле (19) (решение уравнения (18)); 2 - по формуле (33) (решение уравнения (32)); о - численное решение

Анализ результатов расчетов перемещения фронта температурного возмущения д (Ро) показывает, что с увеличением числа приближений время достижения координаты р = 1 уменьшается с

Ро1 = 0,084 (в первом приближении) до Ро* = 0,05 (во втором). И в пределе при п ^да Ро ^ 0, что свидетельствует о приближении к описываемой уравнением (11) бесконечной скорости распространения теплоты. Следовательно, указанное выше противоречие, связанное с определением фронта температурного возмущения в краевой задаче для параболического уравнения теплопроводности, оказывается разрешенным. Необходимость введения фронта возмущения в данном случае связана лишь с упрощением процесса получения приближенного аналитического решения фактически нелинейного дифференциального уравнения с переменными свойствами среды (нелинейность второго рода).

Соотношения (31), (33) представляют решение задачи (11)-(14) во втором приближении. Результаты расчетов по формуле (31) в сравнении с решением задачи (11 )-(14) методом конечных разностей для у = 0,01 приведены на рис. 3. Их анализ показывает, что в диапазоне 0,001 < Ро < 0,05 расхождение результатов не превышает 0,5 %.

Во второй стадии процесса, соответствующей времени Ро < Ро , понятие фронта температурного возмущения теряет смысл. Теплообмен в данном случае происходит по всей толщине пластины, и

+

математическая постановка задачи приводится к виду

5©(р,Ро) _ д_ ЛРо ф (Ро < Ро < да; © (0, Ро) = 1 ; д©(1,Ро)/ Ф = 0.

v(1-р) 5©(р,Ро)

Ф

0 <р< 1);

(34)

(35)

(36)

1,0 ©

0,8 р 1,0

Рис. 3. Распределение температуры в пластине (V = 0,01): Л - по формуле (39) (второе приближение); ° - по формуле (31); --метод конечных разностей

Введем дополнительную искомую функцию

92 (Ро) = ©(1,Ро), (37)

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

Задачи (11)-(14) и (34)-(36) при Ро = Ро1, т. е. при д(Ро) = 1, полностью совпадают. В связи с этим начальное условие для уравнения (34) приводится к соотношению (16), в котором следует положить д (Ро) = 1:

©(р,Ро1) = (1-р)2.

(38)

Решение задачи (34)-(38) принимается в виде

©(р,Ро) = £ ьк Шрк,

(39)

к=0

где Ьк (д2) - неизвестные коэффициенты, определяемые из условий (35)-(37). Подставляя (39), ограничиваясь тремя членами ряда, в (35)-(37) относительно Ь0, Ь, Ь2 будем иметь систему трех алгебраических уравнений, из решения которой находим: Ь0 = 1; Ь1 =-2 (1 - д2); й2 = 1-д2. Соотношение (39) с учетом найденных значе-

ний коэффициентов ЬЛ(д2) (к = 0,1,2) принимает вид

©(р,Ро) = 1 -р(2 -р)(1- 92(ро)) .

(40)

Так как при Ро = Ро д2(Ро) = 0, то соотношение (40) приводится к формуле (38) и, следовательно, уже на этой стадии получения решения начальное условие (38) для функции ©(р,Ро) оказывается выполненным.

Потребуем, чтобы соотношение (40) удовлетворяло не уравнению (34), а некоторому осредненному по пространственной переменной уравнению, т. е. интегралу теплового баланса вида

5©(р, Ро)

ЛРо

dр=í^ Л Лр

dр. (41)

у(1-р) д©(р,Ро)

„ „ 5Р -

Подставляя (40) в (41), относительно неизвестной функции д (Ро) получаем следующее обыкновенное дифференциальное уравнение:

dРo

= -3(д - 1)exp(-v).

(42)

Интегрируя уравнение (42), находим д2 (Ро) = 1 - с, ехр[-3Ро exp(-v)], (43)

где с - константа интегрирования, определяемая из начального условия д2 (Ро) = 0. Формула для нее будет иметь вид С = ехр[3Ро ехр(^)].

Соотношение (40) с учетом (43) принимает вид

©(р,Ро) = 1 - р(2 - р) ехр[-3(Ро - Ро) exp(-v)]. (44)

Соотношение (44) представляет решение задачи (34)-(38) в первом приближении. Оно точно удовлетворяет граничным условиям (35), (36) и начальному условию (38) и приближенно (в первом приближении) - уравнению (34). Отметим, что интеграл теплового баланса (41) решением (44) удовлетворяется точно. Анализ результатов расчетов по формуле (44) в сравнении с расчетом по методу конечных разностей позволяет заключить, что в диапазоне Ро < Ро их расхождение не превышает 6 %.

Для повышения точности решения необходимо увеличивать число членов ряда (39), для определения неизвестных коэффициентов которого будем использовать основные (35)-(37) и дополнительные граничные условия, определяемые таким образом, чтобы искомое решение удовлетворяло уравнению (34) в граничных точ-

0,8

0,6

0,4

0,2

0

0,2

0,4

0,6

ках р = 0 и р = 1. Формула для первого из них, получаемого на основе основного граничного условия (35), совпадает с формулой (22). Еще два дополнительных граничных условия получаются с использованием соотношений (36), (37). Дифференцируя эти соотношения по переменной Ро, получаем:

д 2©(1,Fo)

= 0;

(45)

(46)

дРоф Э0(1,Ро) _ Сд2 (Ро) 5Ро сРо

Продифференцируем уравнение (34) по переменной р и запишем полученное соотношение для точки р = 1:

д2©(1,Ро) _2уд2©(1,Ро) | д3©(1,Ро) (47)

5Роф ф2 ф3 '

Сравнивая (45) и (47), с учетом (36) находим следующее дополнительное граничное условие:

+ = 0. (48)

ф ф

Сравнивая уравнение (34) и соотношение (46), получаем еще одно дополнительное граничное условие в точке р = 1:

д2©(1,Ро) _ Сд2(Ро)

■ (49)

Подставляя (39), ограничиваясь шестью членами ряда в (22), (35)-(37), (48), (49), относительно неизвестных коэффициентов Ьк, (к = 0,5) будем иметь цепочную систему шести алгебраических уравнений: Ь у + 2Ь2 = 0,

(Ь0 + ^р + Ь2р2 + £>3р3 + Ь4р4 + Ь5р5 ) = 1,

Ь + 2Ь + 3Ь + 4Ь + 5Ь = 0, Ь + Ь + Ь + Ь + Ь + Ь = Ъ, (-у)(Ь + 3Ь + 4Ь + 5Ь) = ся2 /СРо, 2Ь у + 3Ь (1 + 2у) +12Ь (1 + у) +10Ь (3 + 2у) = 0^

Ее решение: Ь0 = 1; Ь1 = (д2 Г4 + б0г5 )/ 3/1; Ь2 =у(д2 /4 + б0г5 ^(б/1); Ь = (Я2 (14-у2) + 20 (д2 (1 - у) + у — 2))//1; Ь4 = д (9б + 11у - бу2) +10г4 (д2 - ^/(б/);

b5 = q2(15 + 3v - v2 +12r2(q2 - 1)/(3r1),

где

ri = v-8 ;

2(q2

Г2 = 3-v ;

r3 = 24 - 9v ;

r4 = 9 + 2v ; r5 = q2 - i ; q^ = dq2/cFo.

Подставляя (39) (с учетом найденных значений коэффициентов Ьк, (к = 0,5)) в интеграл теплового баланса (41), относительно неизвестной функции д2 (Ро) будем иметь обыкновенное дифференциальное уравнение

ni

d q? dq9 „„„„ „ 42 1 ~ W2 + 3600 = 0,

dFo2

- + ^2

dFo

(50)

где щ = (66 — v(v + 6))expv ;

ц2 = 1080 + 120v + (540 - 120v) exp v.

Интегрируя уравнение (50), находим q2 (Fo) = i + C exp(|Fo) + C2 exp(|i2Fo), (51) где C, C2 - константы интегрирования;

2 = z + 30(z1 + z2 + z3 + +324)i/2/(v2 exp v + 6(exp v)(11 + v)) ; z = 60v(1 - exp v) + 270exp v + 540 ; Z = 4v(18 + v); z2 = e2v(81- 36v + 4v2) ; z3 = 4ev(15-v2 -3v).

Для определения констант интегрирования используются следующие начальные условия:

q2 (Foi ) = 0 ; dq2 (Foi )/dFo = 0 . (52)

Подставляя (51) в (52), находим Ci = -i/ (i-i2/ ii ); C2 = 12/ (i-i2 ).

Результаты расчетов по формуле (39) во втором приближении (шесть членов ряда) в сравнении с расчетом численным методом при ( v = 0,01) приведены на рис. 3. Из их анализа следует, что в диапазоне 0,025 < Fo <œ расхождение расчетов не превышает 0,5 %. Результаты расчетов по формуле (39) во втором приближении для v = 1,0 в сравнении с решением при v = 0,01 приведены на рис. 4. Их анализ позволяет заключить о существенном влиянии на распределение температуры величины v (расхождение результатов более 25 %). При этом наибольшее расхождение наблюдается в центре симметричной пластины, т. е. там, где расхождение в коэффициенте теплопроводности, определяемом по формуле (5) (где m = v/S ), максимальное. Прикладное значение полученного результата в том, что данный метод можно применять к решению краевых задач с переменными физическими свойствами среды, причем не только к задачам, где коэффициент теплопроводности изменяется монотонно, но и к задачам для многослойных конструкций (композиционные материалы), где он изменяется скачкооб-

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

1,0 о 0,8

чч -О" ЛЧ

\\ \ \ ч \ \ \ \ ^ \ \ ч ч \

\ \ ч N \ \ \ \ \ \ \ ч ч. ч ч - - _ 0,6

\ N \ \ ч ч

ч ч ^0,1

0 0,2 0,4 0,6 0,8 4 1,0

Рис. 4. Распределение температуры в пластине. Расчет по формуле (39) (второе приближение): --V = 0,01; ------- V = 1,0

Результаты расчетов по формуле (31) при V = 0,01 для сверхмалых значений времени (0 < Ро < 10-5) в сравнении с точным решением [1], полученным при v = 0, приведены на рис. 5. Их анализ показывает практическое совпадение результатов двух методов. Этот факт можно объяснить тем, что при сверхмалых значениях времени фронт температурного возмущения проникает на незначительную глубину по координате р (д(Ро) = 2 • 103). В диапазоне столь малой величины пространственной переменной 0<р<2-10"3 происходит столь незначительное изменение коэффициента теплопроводности, что им можно пренебречь (т. е. принять v = 0). Таким образом, уже во втором приближении в первой и второй стадиях процесса получаются достаточно высокой точности приближенные аналитические решения практически во всем диапазоне времени нестационарного процесса.

1,0 0 0,8

0,6

0,4

0,2

5- 10-4

1- 10-3

1,5- 10-3 р 2- 10-3

Рис. 5. Распределение температуры при сверхмалых значениях временной и пространственной переменных (v = 0,01): --по формуле (31); ° - точное решение

Если положить в формулах для (к = 1,2) (являющихся по сути собственными числами краевой задачи (34)-(36)) V = 0 , то получаем следующие их величины: ^ = -2,4709; ц2 = -22,0745 . Точные значения этих чисел: =-2,4674;

ц2 = -22,2066 .

Для нахождения решения задачи (11 )-(14) в третьем приближении необходимо использование еще трех дополнительных граничных условий (4-го, 5-го и 6-го). Для получения четвертого дополнительного граничного условия продифференцируем дополнительное граничное условие (30) по переменной Ро :

а3 ©(р,Ро)

5р25Ро

+ V

а2 ©(р,Ро)

р = 0

араРо

= 0.

(53)

р=0

Дифференцируя уравнение (11), соответственно, один и два раза по переменной р и применяя полученные соотношения для точки р = 0, будем иметь:

а2©(р,Ро)

араРо

а©(р,Ро)

р=0

ар

р=0

-2 V

а2 ©(р,Ро)

ар2

а3 ©(р,Ро)

р=0

а3©(р,Ро)

ар2аРо

ар3

3 а©(р,Ро)

(54)

р=0

р=0

ар

р=0

+3 V2

а2 ©(р,Ро)

ар2

- 3 V

а3 ©(р,Ро)

р=0

ар3

(55)

р=0

а4 ©(р,Ро)

р=0

0

0,6

0,4

0,2

2

= V

+

Подставляя (54) и (55) в (53), получаем четвертое дополнительное граничное условие:

(©V + 4у©'р' + 5у2 ©'р'+ 2у3 ©р) = 0; (56)

где ©^ = д5©/др5 , ©V = д4©/др4;

©''' = д3©/др3; ©|| = д2©/др2; ©'р = д©/др .

Продифференцируем дополнительное граничное условие (33) по Ро:

д3©( р,Ро )

др2дРо

= 0.

(57)

р=д

Дифференцируя уравнение (11) дважды по переменной р и применяя полученные соотношения для точки р = д(Ро), с учетом (28) будем иметь

д3 ©( р,Ро)

др2дРо

= в-д1у

д4 ©( р,Ро)

р=Ч1

д 4

р=д1

-3ув-д1у

д3©( р,Ро)

др 3

(58)

р=д1

Сравнивая соотношения (57) и (58), получаем пятое дополнительное граничное условие:

(©V - 3у©1') = 0. (59)

Для получения шестого дополнительного граничного условия продифференцируем дополнительное граничное условие (28) по переменной Ро:

д 4©( р,Ро)

др3дРо

= 0.

(60)

р=д

Дифференцируя уравнение (11), соответственно, три раза по переменной и применяя полученные соотношения для точки р = д(Ро), с учетом (14), (25) и (28) будем иметь

д4©( р,Ро)

др3дРо

д5©( р,Ро)

р=с/1

д 5

- 4у

р=с/1

д 4©( р,Ро)

(61)

Сравнивая соотношения (60) и (61), получаем шестое дополнительное граничное условие:

(©р - 4у©рУ) = 0 . (62)

Обыкновенное дифференциальное уравнение при у = 0,01 в данном случае будет иметь вид

Сщ 21ехр (-0,01)

д + 375)

С Ро 15б2500д

Ч3 + д2

10б б25 50

+43д +140

(д3 + 2400д2 + 2,28 • 10б д + 5б -107)

+ М + 109 д +1400Л

у10б б25 25 1 у

Из дополнительных граничных условий третьего приближения второй стадии процесса первое условие совпадает с соотношением (56), а второе и третье имеют вид:

©V + 3у©р' + 3у2©Ц = ехр (2у) 4; ©¡V + бу©^ +12у2©I' +10у3©Ц = 0 ,

где = С 2д2/СРо2.

Обыкновенное дифференциальное уравнение относительно д (Ро) при у = 0,01 в третьем приближении будет иметь вид

2,99579д2 (Ро) +1,40159 + 0,07б8905 С д 1

СРо

+0,000882349 - 2,9957 = 0.

сРО2

сРО

Из решения аналогичного уравнения при у = 0 получаются следующие собственные числа: ^ = -2,4б7394 ; ц2 = -22,132366 ; ц3 = -б3,285954. Точное значение третьего собственного числа - ц3 = -б1,б850. Отметим, что в третьем приближении происходит существенное уточнение первого и второго собственных чисел.

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

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

его выполнению и внутри рассматриваемой области.

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

1. Лыков А.В. Теория теплопроводности. - М.: Высш. шк., 1967. - 600 с.

2. Цой П.В. Системные методы расчета краевых задач тепломассопереноса. - М.: Изд-во МЭИ, 2005. - 567 с.

3. Кудинов В.А. Карташов Э.М., Калашников В.В. Аналитические решения задач тепломассопереноса и термоупругости для многослойных конструкций. - М.: Высш. шк., 2005. - 430 с.

4. Канторович Л.В., Крылов В.И. Приближенные методы высшего анализа. - Л.: Физматгиз, 1962. - 708 с.

5. Лыков А.В. Методы решения нелинейных уравнений нестационарной теплопроводности // Энергетика и транспорт. - 1970. - № 5. -С. 109-150.

6. Гудмен Т. Применение интегральных методов в нелинейных задачах нестационарного теплообмена. Проблемы теплообмена: сб. науч. тр. - М.: Атомиздат, 1967. - С. 41-96.

7. Кудинов В.А, Кудинов И.В., Котова Е.В. Дополнительные граничные условия в нестационарных задачах теплопроводности // Теплофизика высоких температур. - 2017. -Т. 55, № 4. - С. 556-563. DOI: 10.1134/S0018151X17040101.

8. Кудинов В.А., Кудинов И.В., Сквор-цова М.П. Обобщенные функции и дополнительные граничные условия в задачах теплопроводности для многослойных тел // Журнал вычислительной математики и математической физики. - 2015. - Т. 55, № 4. - С. 669-680. DOI: 10.1134/S0965542515040089.

9. Кот В.А. Тождества взвешенной температуры // Инженерно-физический журнал. -2015. - Т. 88, № 2. - С. 409-424.

10. Кот В.А. Метод взвешенной температурной функции // Инженерно-физический журнал. - 2014. - Т. 89, № 1. - С. 183-202.

11. Канторович Л.В. Об одном методе приближенного решения дифференциальных уравнений в частных производных // ДАН СССР. - 1934. - Т. 2, № 2. - С. 532-534.

12. Федоров Ф.М. Граничный метод решения прикладных задач математической физики. - Новосибирск: Наука, 2000. - 220 с.

13. Кудинов И. В., Кудинов В.А. Аналитические решения параболических и гиперболических уравнений тепломассопереноса. -М.: Инфра-М, 2013. - 391 с.

14. Кудинов И.В., Кудинов В.А., Котова Е.В. Аналитические решения задач теплопроводности на основе определения фронта теплового возмущения // Изв. вузов. Математика. - 2016. - № 11. - С. 27-41.

References

1. Lykov, A.V. Teoriya teploprovodnosti [Theory of heat conductivity]. Moscow: Vysshaya shkola, 1967. 600 p.

2. Tsoy, P.V. Sistemnye metody rascheta kraevykh zadach teplomassoperenosa [System methods of calculation of boundary value problems of heat and mass transfer]. Moscow: Iz-datel'stvo MEI, 2005. 567 p.

3. Kudinov, V.A. Kartashov, E.M., Kalashni-kov, V.V. Analiticheskie resheniya zadach tep-lomassoperenosa i termouprugosti dlya mnog-osloynykh konstruktsiy [Analytical solutions to problems of heat and mass transfer and thermoe-lasticity for multilayer structures]. Moscow: Vysshaya shkola, 2005. 430 p.

4. Kantorovich, L.V., Krylov, V.I. Priblizhen-nye metody vysshego analiza [Approximate methods of higher analysis]. Leningrad: Fizmatgiz, 1962. 708 p.

5. Lykov, A.V. Metody resheniya nelineynykh uravneniy nestatsionarnoy teploprovodnosti [Methods for solving non-linear equations of transient heat conduction]. Energeti-ka i transport, 1970, no. 5, pp. 109-150.

6. Gudmen, T. Primenenie integral'nykh metodov v nelineynykh zadachakh nestatsionar-nogo teploobmena [Application of integral methods in nonlinear problems of unsteady heat transfer]. Sbornik nauchnykh trudov «Problemy teploo-bmena» [Collection of scientific works «Problems of heat transfer»]. Moscow: Atomizdat, 1967, pp. 41-96.

7. Kudinov, V.A, Kudinov, I.V., Kotova, E.V. Dopolnitel'nye granichnye usloviya v nes-tatsionarnykh zadachakh teploprovodnosti [Additional boundary conditions in transient heat conduction problems]. Teplofizika vysokikh temperature, 2017, vol. 55, no. 4, pp. 556-563. DOI: 10.1134/S0018151X17040101.

8. Kudinov, V.A., Kudinov, I.V., Skvort-sova, M.P. Obobshchennye funktsii i dopolnitel'nye granichnye usloviya v zadachakh teploprovodnosti dlya mnogosloynykh tel [Generalized functions and additional boundary conditions in heat-conductivity problems for multilayer bodies]. Zhur-nal vychislitel'noy matematiki i matematicheskoy fiziki, 2015, vol. 55, no. 4, pp. 669-680. DOI: 10.1134/S0965542515040089.

9. Kot, V.A. Tozhdestva vzveshennoy tem-peratury [Identities of weighted temperature]. In-zhenerno-fizicheskiy zhurnal, 2015, vol. 88, no. 2, pp. 409-424.

10. Kot, V.A. Metod vzveshennoy tempera-turnoy funktsii [Method of weighted temperature function]. Inzhenerno-fizicheskiy zhurnal, 2014, vol. 89, no. 1, pp. 183-202.

11. Kantorovich, L.V. Ob odnom metode priblizhennogo resheniya differentsial'nykh uravneniy v chastnykh proizvodnykh [On one

method of approximate solution to partial differential equations]. DAN SSSR, 1934, vol. 2, no. 2, pp. 532-534.

12. Fedorov, F.M. Granichnyy metod resheniya prikladnykh zadach matematicheskoy fiziki [A boundary method of solving applied problems of mathematical physics]. Novosibirsk: Nauka, 2000. 220 p.

13. Kudinov, I. V., Kudinov, V.A. Analitich-eskie resheniya parabolicheskikh i giperbolich-eskikh uravneniy teplomassoperenosa [Analytical

solutions to parabolic and hyperbolic equations of heat and mass transfer]. Moscow: Infra-M, 2013. 391 p.

14. Kudinov, I.V., Kudinov, V.A., Kotova, E.V. Analiticheskie resheniya zadach teploprovodnosti na osnove opredeleniya fronta teplovogo vozmushcheniya [Analytic solutions to heat transfer problems based on determining the temperature perturbation front]. Izvestiya vuzov. Matemat-ika, 2016, no. 11, pp. 27-41.

Котова Евгения Валериевна,

ФГБОУВО «Самарский государственный технический университет», кандидат технических наук, доцент кафедры теоретических основ теплотехники и гидромеханики, телефон (846) 332-42-35, e-mail: [email protected] Kotova Yevgenia Valerievna,

Samara State Technical University, Candidate of Engineering Sciences (PhD), Associate Professor of the Department of Theoretical Foundations of Heat Engineering and Hydromechanics, telephone (846) 332-42-35, e-mail: [email protected]

Еремин Антон Владимирович,

ФГБОУВО «Самарский государственный технический университет», кандидат технических наук, доцент кафедры теоретических основ теплотехники и гидромеханики, e-mail: [email protected] Eremin Anton Vladimirovich,

Samara State Technical University, Candidate of Engineering Sciences (PhD), Associate Professor of the Department of Theoretical Foundations of Heat Engineering and Hydromechanics, e-mail: [email protected]

Кудинов Василий Александрович,

ФГБОУВО «Самарский государственный технический университет», доктор физико-математических наук, профессор, зав. кафедрой теоретических основ теплотехники и гидромеханики, телефон (846) 332-42-35, e-mail: [email protected] Kudinov Vasily Aleksandrovich,

Samara State Technical University,Doctor of Physical and Mathematical Sciences (Post-doctoral degree), Professor, Head of the Department of Theoretical Foundations of Heat Engineering and Hydromechanics, telephone (84б) 332-42-35, e-mail: [email protected]

Ткачев Василий Константинович,

ФГБОУВО «Самарский государственный технический университет», аспирант кафедры теоретических основ теплотехники и гидромеханики, телефон (846) 332-42-35, e-mail: [email protected] Tkachev Vasily Konstantinovich,

Samara State Technical University, Post-graduate student of the Department of Theoretical Foundations of Heat Engineering and Hydromechanics, telephone (846) 332-42-35, e-mail: [email protected]

Кузнецова Анастасия Эдуардовна,

ФГБОУВО «Самарский государственный технический университет», кандидат технических наук, старший преподаватель кафедры теоретических основ теплотехники и гидромеханики, телефон (846) 332-42-35, e-mail: [email protected] Kuznetsova Anastasia Eduardovna,

Samara State Technical University, Candidate of Engineering Sciences (PhD), Senior Lecturer of the Department of Theoretical Foundations of Heat Engineering and Hydromechanics, telephone (846) 332-42-35, e-mail: [email protected]