Получение приближенных аналитических решений гиперболических уравнений теплопроводности

443100, Россия, г. Самара, ул. Молодогвардейская, 244

E-mail: totig@yandex.ru

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

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

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

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

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

(1)

где  - температура;

 - координата;

 - время;

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

 - температура стенки;

 - толщина пластины;

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

 - коэффициент релаксации (время релаксации);

 - скорость продвижения тепловой волны.

Введем следующие обозначения:

; ;

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

 - безразмерная координата;

 - число Фурье;

- безразмерный коэффициент релаксации.

С учетом принятых обозначений задача (1) приводится к виду

; (2)

; (3)

; (4)

; (5)

, (6)

где

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

, (7)

()

где ();

(); ;

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

Исследования решения (7) показали, что для некоторых чисел после достижения фронтом тепловой волны центра пластины наблюдается обратная тепловая волна, также имеющая скачок температуры на ее фронте [7, 9].

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

. (8)

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

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

. (9)

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

. (10)

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

Следует отметить, что сходимость ряда (7) существенно зависит от величины числа Фурье, для которого необходимо определять распределение температуры по толщине пластины. Например, в диапазоне числа Фурье для сходимости достаточно всего нескольких членов этого ряда. При сходимость имеет место при членах ряда. Для всех чисел Фурье, при которых происходит скачок температурных кривых, число членов ряда (7), необходимых для его сходимости, существенно возрастает

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

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

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

; (11)

; (12)

, (13)

где . (14)

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

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

, (15)

где  - неизвестная функция времени;

 - неизвестный коэффициент.

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

. (16)

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

, (17)

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

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

(18)

где ; ; ; ;; ;

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

 - это время, в течение которого температура на фронте тепловой волны становится равной единице.

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

(19)

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

, (20)

где .

После определения коэффициента k из (20) решение задачи (2), (11) - (14) находится из (17). Анализ результатов расчетов по формуле (17) при в сравнении с точным решением, определяемым по формуле (7), позволяет заключить, что их отличие на участке временной координаты не превышает 2,5 %, причем с уменьшением числа точность решения по формуле (17) возрастает и уже при оно практически совпадает с точным. Величина коэффициента для оказалась равной k = 2,5376. Величина , найденная из точного аналитического решения (7), составляет .

Результаты расчетов по формуле (17) для даны на рис. 1. Их анализ приводит к выводу о том, что отклонение от точного решения в данном случае составляет около 1,5 %. Следовательно, с уменьшением величины точность решения по формуле (17) возрастает. Значения величин k и для составляют: ; .

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

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

Выражая в решении (17) координату как функцию и , получаем следующую формулу для построения изотерм (при ):

. (21)

Задаваясь произвольными значениями , из (21) получаем графики изотерм в координатах . Графики изотерм, полученные по формуле (21), в сравнении с кривыми изотерм, полученными из точного аналитического решения, даны на рис. 2. Штриховой линией обозначена линия перемещения фронта тепловой волны по координате во времени . Отметим, что по формуле (21) построены изотермы лишь до момента времени . Анализ движения изотерм по пространственной координате во времени показывает, что каждая из них возникает внутри тела в строго определенный момент времени и в строго определенной точке по координате . Отличие изотерм, полученных по формуле (21), от изотерм, найденных на основе использования точного решения, находится в пределах 1,5 %.

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

. (22)

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

Определяя первые производные по времени от соотношений (21), (22), можно найти безразмерные скорости движения изотерм по координате во времени. Например, при формула для определения скоростей изотерм имеет вид

уравнение гиперболический теплопроводность

, (23)

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

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

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

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

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

, (24)

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

Записывая соотношение (24) для точек 1, 2, 3, 4 кривой, считая, что температуры в этих точках наблюдаются соответственно для чисел Фурье, равных ; ; ; , для определения неизвестных коэффициентов будем иметь систему четырех алгебраических линейных уравнений, в которой значения полученных из решения температур в точках 1, 2, 3, 4, были соответственно равны ; ; ; . Из решения этой системы уравнений находим

; ; ; (25)

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

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

(26)

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

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

Рис. 1. Изменение температуры в пластине ( ). - по формуле (7) (точное решение); - - - - - - по формуле (17)

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

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

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

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

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

Рис. 3. Изменение скоростей движения изотерм по координате во времени : - - - - - - линия перемещения фронта тепловой волны во времени

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

- точное решение; - по формуле (24)

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