Получение приближенных аналитических решений гиперболических уравнений теплопроводности Текст научной статьи по специальности «Физика»

ВЕСТН. САМАР. ГОС. ТЕХН. УН-ТА. СЕР. ТЕХНИЧЕСКИЕ НАУКИ. 2017. № 4 (56)

УДК 536.2 (075)

ПОЛУЧЕНИЕ ПРИБЛИЖЕННЫХ АНАЛИТИЧЕСКИХ РЕШЕНИЙ ГИПЕРБОЛИЧЕСКИХ УРАВНЕНИЙ ТЕПЛОПРОВОДНОСТИ

И.В. Кудинов, О.Ю. Курганова, Г.Н. Максименко, Л.С. Абишева

Самарский государственный технический университет 443100, Россия, г. Самара, ул. Молодогвардейская, 244 E-mail: [email protected]

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

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

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

Игорь Васильевич Кудинов (к.т.н., доц.), доцент кафедры «Теоретические основы теплотехники и гидромеханики».

Ольга Юрьевна Курганова, аспирант.

Любовь Сергеевна Абишева, ассистент кафедры «Теоретические основы теплотехники и гидромеханики».

Галина Николаевна Максименко, аспирант.

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

Математическая постановка задачи теплопроводности для бесконечной пластины при симметричных граничных условиях первого рода с учетом конечной скорости распространения теплового возмущения имеет вид [1-3, 6-9]

3Fo ■""r 3Fo2 ЭХ2

в(Х,0)=1. ?

Эв (Х,0)/ ЭFo = 0.

(1)

д ?(х, т) д 2t (х, т) д 2t (х, т) , . .

; + тг—^^=а—(т>0; о<х<8)

дт дт2 дх2

, т д t (х ,0) п д t(0, т) п . t(х,0) = tо; ' = 0; —^^ = 0; t(8,т) = tcг,

дт дх

где t - температура; х - координата; т - время;

t0 - начальная температура;

- температура стенки; 5 - толщина пластины;

а - коэффициент температуропроводности; 2

тг = а / м> - коэффициент релаксации (время релаксации); ^ - скорость продвижения тепловой волны. Введем следующие обозначения:

© = С - tсг) /& - tсг); Х = х/ 5 ; Р° = а т 7 55; = а 52 ,

где © - относительная избыточная температура; X - безразмерная координата; Бо - число Фурье;

Б°г - безразмерный коэффициент релаксации. С учетом принятых обозначений задача (1) приводится к виду

дв(^о) + р° д 2в(Х,Р°) = д 2е(Х,Б°) . (Р° > 0; 0 <Х< 1) (2)

(3)

(4)

Эв (0,Fo)/ ЭХ = 0. (5)

(6)

в (1,Fo) = 0

где For = a t = const .

В работах [7, 9] путем совместного использования метода разделения переменных и ортогонального метода Л.В. Канторовича получено точное аналитическое решение задачи (2) - (6), которое имеет вид

0(X,Fo) = £ [Clk exp(zuFo) + C2k exp(Z2kFo)]cos I rPX I, (7)

k=l V 2 J

(r = 2k -1; k = ) где zlk = (-1 + 71 - 4Forvk)/(2For) (i = 1,2; k = ); vк = r2p2/4

4(-1)k z2k >

(r = 2k -1; k =1 ¥); C1k =-C2kz2k/z1k; C2k =-/

rp

1 —

V z1k J

Результаты расчетов безразмерной температуры по формуле (7) Fo r = 6,25 • 10 -3 даны на рис. 1. Их анализ позволяет заключить, что при малых значениях безразмерного времени Fo на фронте температурного возмущения наблюдается скачок температурных кривых, то есть по существу образуется фронт тепловой волны, на границе которого наблюдается скачок температуры от ее значения в точке скачка до величины начальной температуры. Следовательно, область, находящаяся за пределами фронта тепловой волны, оказывается невозмущенной и имеющей начальную температуру.

Исследования решения (7) показали, что для некоторых чисел For 2

(For >10"2) после достижения фронтом тепловой волны центра пластины наблюдается обратная тепловая волна, также имеющая скачок температуры на ее фронте [7, 9].

Анализ полученных результатов позволил заключить, что температура в точках скачка (на фронте тепловой волны) описывается формулой

t(t) = t0 -exp[-t/(2tr)]. (8)

Эта формула полностью согласуется с соотношением для температуры на фронте тепловой волны, полученным в работах [1, 3, 4] для полупространства.

Соотношение (8) в безразмерном виде будет

0(Fo) = 1 - exp (- 0,5FoFor1) (9)

Расчеты, выполненные по формуле (7), показали линейную закономерность движения фронта тепловой волны по координате х во времени Хф = 5 - wt , или

в безразмерном виде

Хф (Fo) = 1 - FoFo-05. (10)

Анализ результатов расчетов по формуле (7) позволяет заключить, что после того как на фронте тепловой волны температура становится равной начальной (для For = 6,2540 3 это происходит при Fo = 0,073), для всех последующих моментов времени вплоть до достижения стационарного состояния распределение температуры, найденное по формуле (7), полностью совпадает с решением параболического уравнения теплопроводности при тех же граничных условиях.

Следует отметить, что сходимость ряда (7) существенно зависит от величины числа Фурье, для которого необходимо определять распределение температуры по толщине пластины. Например, в диапазоне числа Фурье 0,1 < Fo < ¥ для сходимости достаточно всего нескольких членов этого ряда. При 0,0001 < Fo < 0,1 сходимость имеет место при 10 ^1000 членах ряда. Для всех чисел Фурье, при которых происходит скачок температурных кривых, число членов ряда (7), необходимых для его сходимости, существенно возрастает

(от 104 при Бо = 10-5 до 106 при Бо=10- ). При дальнейшем уменьшении числа Фурье количество членов ряда (7) может достигать нескольких миллионов.

Учитывая указанные выше трудности получения аналитических решений для сверхмалых значений времени, следует отметить, что актуальной является проблема нахождения более простых выражений, описывающих температурное состояние для моментов времени, при которых происходит скачок температурных кривых. С этой целью рассмотрим задачу (2) - (6) лишь для первой стадии процесса 0 < Бо < Б°1, где Б°1 - время, при котором температура на фронте тепловой волны становится равной начальной температуре, т. е. ©(Хф ^о^ = 1.

Математическая постановка задачи в данном случае включает уравнение (2) при (0 < Бо < Бо1; 0 < X < Хф (Бо)), а также следующие начальные и граничные условия:

© а;°) =0; (11)

®(1,Ро) = 0; (12)

в(Хф ,Бо) = 1 - ехр (- 0,5РоРо-1)

(13)

где ХФ(Ро) = 1 - РоРо -0". (14)

Отметим, что в задаче (2), (11) - (14) учтены соотношения (9), (10) и к тому же в отличие от (3) принимается нулевое начальное условие. Это связано с тем, что согласно соотношению (13), совпадающему с соотношением (9), температура на фронте тепловой волны принимается известной в любой момент времени первой стадии процесса. Таким образом, начальное условие (11) следует из соотношения (13), согласно которому при Бо=0 ©(1;0) = 0 .

Следуя методу Л.В. Канторовича, решение задачи (2), (11) - (14) принимаем в виде произведения двух функций

©(Х,Ро) = Ь(?о)(1 -X*), (15)

где Ь (Бо) - неизвестная функция времени; * - неизвестный коэффициент. Соотношение (15) удовлетворяет граничному условию (12). Для нахождения неизвестной функции Ь (Бо) используем граничное условие (13). Подставляя (15) в (13), получаем

Ь(Бо) = [1 - ехр (-0,5ЕоБо-1)] /(1 - £|). (16)

Подставляя (16) в (15), будем иметь

1 -Х*п пс^-ь

:к

ф (17)

©(Х,Бо) = —X [1 - ехр(-0,5РоРо-1)]

1 ХФ

где Хф находится из (14).

Соотношение (17) независимо от величины коэффициента к удовлетворяет начальному условию (11) и граничным условиям (12), (13). Неизвестный коэффициент к находится так, чтобы удовлетворялось уравнение (2). Подставляя (17) в (2), относительно коэффициента к получаем следующее трансцендентное уравнение:

0,5A7 _ kA1 + 2AJc2 _ kA3 + Atk2 _ kA, _ 0,25A7 _ Q

A3 F°r A4A2 A7A2( } A5A2 A6A2 A6A2 F°rA3 (¡g)

где 4_1-A3; A2 _ 1-Fo/Fc^5; A3 _exp(_0,5Fo/For); A4 _A2F°0'5;

A5 _ A2F°,'5; As _ AfFo, ; A7 _ 1_ A.

Анализ уравнения (18) приводит к заключению об отсутствии в нем переменной X, что связано с выполнением соотношений (13), (14), полученных из точного аналитического решения вида (7) задачи (2) - (6). Таким образом, в уравнении (18) коэффициент k зависит лишь от времени Fo. Для определения его числового значения проинтегрируем уравнение (18) по переменной Fo в пределах от Fo _ 0 до Fo _ Foj, где величина Foj для каждого конкретного значения Fo r находится из соотношения (7), положив 0(X,Fo1) _ 1, то есть Fo1 - это время, в течение которого температура на фронте тепловой волны становится равной единице.

Определяя интеграл невязки уравнения ( 18), получаем

0,5 A7 kA 2 Ak2 kA3 Ark (k _ 1) 0,25 A7

-+-~—T~--Ц- +-

4_1For A4 A2_k A7Af_k) A5 A2_k AsA2_k Fo2 A3"1

dFo = 0.

(19)

Определяя интегралы в (20), относительно коэффициента к получаем следующее соотношение:

к^о, ,Бо1) = 2Мо0 '^/^(Я - 2) -

B _ ln[(For _ 2Fo0'5Fo1 + Fof)/For].

(20)

где B _ in[(F or _ 2F o.• ~п~2л

После определения коэффициента k из (20) решение задачи (2), (11) - (14) находится из (17). Анализ результатов расчетов по формуле (17) при

For _ 6,25 10_3 в сравнении с точным решением, определяемым по формуле (7), позволяет заключить, что их отличие на участке временной координаты 0 < Fo < Foj не превышает 2,5 %, причем с уменьшением числа Fo точность решения по формуле (17) возрастает и уже при Fo < 0,02 оно практически совпадает с точным. Величина коэффициента k для For _ 6,25 10_3 оказалась равной k = 2,5376. Величина Fob найденная из точного аналитического решения (7), составляет Fq _ 0,073.

_7

Результаты расчетов по формуле (17) для For _ 10 даны на рис. 1. Их анализ приводит к выводу о том, что отклонение от точного решения в данном случае составляет около 1,5 %. Следовательно, с уменьшением величины For точность решения по формуле (17) возрастает. Значения величин k и Foj для For _ 10_7 составляют: k _ 531,8699 ; Fo1 _ 9,48683298 ■ 10_6 .

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

теплопроводности идентифицирован коэффициент релаксации тг , практическое определение которого представляет значительные трудности [2, 3].

Ввиду достаточно простой зависимости безразмерной температуры © от пространственной переменной X соотношение (17) удобно использовать для построения изотерм в координатах X - Бо . Отметим, что построение изотерм путем непосредственного выражения координаты Х через переменную Бо из решения (7) не представляется возможным. В связи с этим изотермы точного решения строились путем многовариантных расчетов по формуле (7) с последующим графическим построением кривых для каждой отдельной изотермы.

Выражая в решении (17) координату X как функцию © (X ,Бо) и Бо, получаем следующую формулу для построения изотерм (при Бо г = 6,25 ■ 10-3):

2 537586 "10,394075

1 - (1 - 12,649110Ро)

X(0,Fo)=

1 -0(X,Fo)-

(21)

1 - exp(-80Fo)

Задаваясь произвольными значениями 0 = const, из (21) получаем графики изотерм в координатах X-Fo . Графики изотерм, полученные по формуле (21), в сравнении с кривыми изотерм, полученными из точного аналитического решения, даны на рис. 2. Штриховой линией обозначена линия перемещения фронта тепловой волны по координате X во времени Fo. Отметим, что по формуле (21)

построены изотермы лишь до момента времени Fo = Fc^. Анализ движения изотерм по пространственной координате во времени показывает, что каждая из них возникает внутри тела в строго определенный момент времени и в строго определенной точке по координате Х . Отличие изотерм, полученных по формуле (21), от изотерм, найденных на основе использования точного решения, находится в пределах 1,5 %.

Формула для построения изотерм при малых значениях For (For =10-7) приводится к виду

-.0,00188016

X(0,Fo) =

1 -0(X,Fo)

1 - (1 - 3162,27766Fo)

531,8699

(22)

1 - ехр(-5000000Бо)

Анализ графиков изотерм, найденных по формуле (22), позволяет заключить, что при малых значениях Бог в диапазоне чисел 0 < Бо < Бо1 изотермы принимают вид практически прямых линий, слабо наклоненных в горизонтальном направлении. Следовательно, скорости движения изотерм в этом интервале времени будут незначительными.

Определяя первые производные по времени Бо от соотношений (21), (22), можно найти безразмерные скорости движения изотерм по координате X во вре-

мени. Например, при For = 6,25 10 терм имеет вид

dX _ 0,39075

-3

формула для определения скоростей изо-

u(0,Fo)=

dFo X(Q, Fo)'

0,6059

- 32,0982040

(1 - 12,64911Fo) 1 - exp(-80Fo)

1,537586

- +

+ 800

1 - (1 -12,64911Fo)

2,537586

1 - exp(-80Fo)

где X(Q,Fo) определяется по формуле (21).

-exp(-80Fo)

(23)

Результаты расчетов скоростей движения изотерм по формуле (23) даны на рис. 3 (скорости изотерм для Fo > 0,073 получены с использованием точного решения (7)). Отрицательные знаки скоростей объясняются тем, что изотермы движутся противоположно направлению оси X . Анализ полученных результатов позволяет заключить, что изотермы, возникая внутри тела в определенные моменты времени в определенных точках координаты X, уже в момент возникновения имеют бесконечно большие начальные скорости, что объясняется бесконечно большой величиной теплового потока на фронте тепловой волны. С увеличением времени скорости уменьшаются до некоторого минимального для каждой изотермы значения. При дальнейшем увеличении времени скорости изотерм начинают возрастать (по абсолютной величине), устремляясь к бесконечным значениям при приближении изотерм к граничной (адиабатной) стенке (X = 0), где в исходной задаче задано граничное условие (5). Отметим, что применительно к рис. 3 для любой изотермы 0 = const координата X , для которой наблюдается соответствующая скорость движения изотерм, находится из графиков рис. 2.

Преимущества полученного выше аналитического решения вида (17) состоят в том, что при использовании экспериментальных данных по изменению во времени температуры в какой-либо одной точке пластины путем решения обратной задачи теплопроводности может быть найдена величина числа For = атг/8г, из

которого можно найти коэффициент релаксации tr.

Допустим, что из эксперимента известно изменение температуры в точке X = 0,8 в диапазоне числа Фурье Fo1 < Fo < Fo4 (в качестве экспериментальных данных будем использовать результаты точного решения задачи (2) - (6) вида (7)). Кривая изменения полученной таким путем температуры дана на рис. 4.

Аппроксимируем эту кривую следующей функцией:

3

0(0,8;Fo) = £btFoi, (24)

i=0

где bi (i = 0,3)- неизвестные коэффициенты, для определения которых используются расчетные данные, полученные по формуле (7).

Записывая соотношение (24) для точек 1, 2, 3, 4 кривой, считая, что температуры в этих точках наблюдаются соответственно для чисел Фурье, равных Fo1 = 0,005 ; Fo2 = 0,02 ; Fo3 = 0,04 ; Fo4 = 0,06 , для определения неизвестных

коэффициентов bi (i = 0,3) будем иметь систему четырех алгебраических линейных уравнений, в которой значения полученных из решения температур в точках 1, 2, 3, 4, были соответственно равны 0(Fo1) = 0,933; 0(Fo2) = 0,659; 0(Fo3) = 0,497 ; Q(Fo4) = 0,441. Из решения этой системы уравнений находим

b1 = 1,064 . b2 = -28,593 . b3 = 475,705 . b4 = -2873,298. (25)

Результаты расчетов по формуле (24) с учетом (25) приведены на рис. 4. Их анализ позволяет сделать заключение о практическом совпадении результатов аппроксимации с точным аналитическим решением, имитирующим экспериментальные данные.

Из решения обратной задачи с использованием соотношения (17) можно идентифицировать (восстановить) число Fo r . Подставляя (24) в левую часть ре-

шения (17), положив X, = 0,8 и определяя интеграл от полученного соотношения в пределах Бо £ Бо£ Бод, получаем

Г (1,064 - 28,593Бо + 475,705Бо2 - 2873,298Бо3 ^Бо =

Бо,

Ро4 1 - о 8^ 1

= Г-^г-1 - ехр^^оБо^Ш Бо,

1 - БоБо0,5 :

Бо^ :

(26)

где к находится из соотношения (20).

Интеграл в правой части может быть найден лишь путем численных расчетов, в результате которых получаем Бо: = 0,00624. Точная величина Бо:, при которой было получено аналитическое решение вида (7), Бо: = 0,00625 . Таким образом, погрешность аппроксимации составляет 0,002 %. В случае использования экспериментальных данных при решении обратных задач теплопроводности точность идентификации будет зависеть от точности выполнения эксперимента.

1,0 0 0,9

0,8

0,7

0,6

0,5

0,4

0,3

0,2

0,1

NN. "V.

О- чч Ч\

Ч^-

V Ч Ч хь Ч 4 ч\

ч. ч NN X XX ч^ \\ Ч Ч ч-

V ч\Л \\\\ Чх ч^

Чх

\

0,9984

0,9986

0,9988

0,999

0,9992

0,9994

0,9996

0,9998

X

1,0

-7 -

Рис. 1. Изменение температуры в пластине (Бо: = 10 ; = 9,48683298 10" к = 531,8699). - - по формуле (7) (точное решение);------по формуле (17)

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

>-6 .

1,0 0,9 0,8 0,7 0,6 0,5 0,4 0,3 0,2 0,1

0 0,04 Foi 0,12 0,16 0,2 0,24 0,28 0,32 0,36 Fo 0,4

Рис. 2. Распределение изотерм 0 = const по координате £ во времени:

.....- линия перемещения фронта тепловой волны.

- - по формуле (7);

о - по формуле (21); For = 6,25 10-3; Foj = 0,073

Рис. 3. Изменение скоростей движения изотерм по координате X во времени Бо:

Fo = 6,25 10

-3 .

■ - линия перемещения фронта тепловой волны во времени

о

4

'О,. 3

__4 С

0,01

0,02

0,03

0,04

0,05 Fo 0,06

Рис. 4. Результаты аппроксимации точного решения соотношением (24): --точное решение; 0 - по формуле (24)

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Шашков А.Г., Бубнов В.А., Яновский С.Ю. Волновые явления теплопроводности: системно-структурный подход. Изд. 2-е, доп. - М.: Едиториал УРСС, 2004. - 296 с.

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

3. Лыков А.В. Тепломассоперенос. - М.: Энергия, 1971. - 560 с.

4. Баумейстер К., Хамилл Т. Гиперболическое уравнение теплопроводности. Решение задачи о полубесконечном теле // Теплопередача. - 1969. - № 4. - С. 112-119.

5. ЖоуД., Касас-БаскесХ., ЛебонДж. Расширенная необратимая термодинамика. - М.Ижевск: НИЦ «Регулярная и хаотическая динамика»: Институт компьютерных исследований, 2006. - 528 с.

6. Карташов Э.М. Аналитические методы в теории теплопроводности твердых тел. - М.: Высшая школа, 2001. - 550 с.

7. Кудинов В.А., Кудинов И.В. Об одном методе получения точного аналитического решения гиперболического уравнения теплопроводности на основе использования ортогональных методов // Вестник Самарского государственного технического университета. Сер. Физико-математические науки. - 2010. - № 5 (21). - С. 159-170.

8. Кудинов В.А., Кудинов И.В. Получение и анализ точного аналитического решения гиперболического уравнения теплопроводности для плоской стенки // Теплофизика высоких температур. - 2012. -Т. 50. - № 1. - С. 118-125.

9. Кудинов В.А., Кудинов И.В. Исследование теплопроводности с учетом конечной скорости распространения теплоты // Теплофизика высоких температур. - 2013. - Т. 51. - № 2. -С. 301-310.

Статья поступила в редакцию 3 октября 2017 г.

OBTAIN APPROXIMATE ANALYTICAL SOLUTIONS OF HYPERBOLIC EQUATIONS THERMAL CONDUCTIVITY

I. V. Kudinov, O. Y. Kurganova, G.N. Maksimenko, L.S. Abisheva

Samara State Technical University

244, Molodogvardeiskaya st., Samara, 443100, Russian Federation E-mail: [email protected]

Based on the use of orthogonal method L. V. Kantorovich obtained an approximate analytical solution of the hyperbolic heat equation for an infinite plate with symmetric boundary conditions of the first kind. Found solution describes the thermal state of a plate on the first stage of the process, i.e. before the time when the temperature at the front of the shock of a thermal wave is equal to the initial temperature. The solution is very simple works algebraic coordinate functions to an exponential function of time, which allows the study of the thermal state of the body in the fields of isotherms with the determination of the velocity of their motion on spatial coordinate in time. The official decision on the experimental values of the temperature in any point of the plate by solving the inverse task allows you to find the coefficient of relaxation, experimental determination of which is difficult.

Keywords: hyperbolic equation, finite speed of heat propagation, analytical solution, isotherms, speeds of movement of isotherms, relaxation coefficient, inverse heat conduction problem.

Igor V. Kudinov (Ph.D. (Techn.)), Associate Professor. Ol 'ga Y. Kurganova, Postgraduate Student. Lubov' S. Abisheva, Postgraduate Student. Galina N. Maksimenko, Postgraduate Student.