Введение
Перспективным методом экологичного и экономичного получения водорода является пиролиз (термическое разложение без доступа кислорода) метана в жидкометаллических средах (олово, висмут и др.). При его практическом применении важной проблемой является появление в керамическом реакторе, содержащем жидкий металл, дефектов в виде трещин, нарушающих его герметичность. Трещинообразование связано с температурными напряжениями, возникающими вследствие больших градиентов температур в процессе пусков и остановов реактора. Для выяснения причин их появления необходимо выполнить исследование температурного состояния реактора в указанных процессах. Решению в данном случае подлежит задача теплопроводности при несимметричных граничных условиях третьего рода. Ее точное аналитическое решение имеет вид бесконечного ряда, собственные числа которого находятся из трансцендентных уравнений, решаемых лишь численно. Получаемое в этом случае решение исходной задачи затруднительно применять для инженерных приложений. В настоящей работе для решения таких задач используется приближенный аналитический метод, практически эквивалентный точному, основанный на использовании дополнительных граничных условий. Основной проблемой, ограничивающей использование приближенных методов (Л. В. Канторовича, Галеркина, интегрального метода и др.), является недостаточная точность решений из-за ограничения числа приближений по причине плохой обусловленности матриц систем алгебраических уравнений [1, с. 101–117; 2, с. 114–125; 3]. Подобные системы возникают при определении неизвестных коэффициентов аналитических решений, принимаемых в форме алгебраических или тригонометрических полиномов с целью выполнения уравнения и краевых условий задачи. Универсальность приближенных методов состоит в выполнении не исходного дифференциального уравнения, а его осредненного аналога – интеграла теплового баланса, в связи с чем эти методы могут быть применены к задачам, решения которых классическими точными аналитическими методами (Фурье, интегральных преобразований и др.) не могут быть получены [4, с. 125–132; 5; 6].
В ряде работ [3; 4, с. 125–132; 5] для повышения точности приближенных методов используются дополнительные граничные условия (ДГУ), оказывающие определяющее влияние на обусловленность матриц, что связано со следующими обстоятельствами. Дополнительные граничные условия приводят к выполнению уравнения в граничных точках [7–9]. В работах [4, с. 125–132; 5] (в результате доказательства теорем) утверждается, что выполнение уравнения на границах при большом числе ДГУ приводит к его выполнению и внутри области.
Используемые в указанных работах ДГУ выполняются при любом другом методе получения решения [10; 11], в том числе и в точных аналитических методах (в чем можно убедиться подстановкой соответствующих решений в ДГУ), однако в этих методах они не выделяются в качестве отдельного рассмотрения [12; 13]. Эффективность их использования, как будет показано ниже, связана с тем, что ДГУ представляются в виде производных высокого порядка от уравнения задачи Штурма – Лиувилля. Удовлетворение таких ДГУ в граничной точке ξ = 0 приводит к цепочным системам алгебраических линейных уравнений относительно неизвестных коэффициентов решений, в которых каждое уравнение содержит лишь одно неизвестное. Такие системы имеют сильно разреженные матрицы коэффициентов ленточного типа с минимальной шириной ленты, проблемы плохой обусловленности которых вообще не возникают. В отличие от [4, с. 125–132; 5], в настоящей работе доказательство теоремы связано с решением не исходного уравнения в частных производных, а уравнения краевой задачи Штурма – Лиувилля, дополнительные граничные условия для которого выполняются лишь в одной граничной точке ξ = 0.
Цель исследования – разработка метода получения высокоточного решения задачи теплопроводности для пластины с неоднородными условиями третьего рода путем использования ДГУ, приводящих к граничному выполнению уравнения, что обеспечивает его выполнение и внутри области, исключая интегрирование по пространственной координате.
Материалы и методы исследования
Найдем решение задачи теплопроводности для керамической стенки реактора пиролиза метана с несимметричными условиями третьего рода. Снаружи реактор омывается горячими газами печи – греющая среда, а внутри через слой расплавленного металла, находящего в начальном диапазоне времени в гранулированном состоянии, движется нагреваемый метан. Математическая постановка задачи:
; (τ > 0; 0 < x < δ); (1)
; (2)
; (3)
, (4)
где t, x, τ– температура, координата, время; t0 – начальная температура; t1 – температура горячих газов; t2 – температура метана; α1, α2 – коэффициенты теплоотдачи; λ – коэффициент теплопроводности; a – коэффициент температуропроводности; δ – толщина стенки реактора.
Обозначим:
;
;
;
;
,
где T, ξ, Fo – безразмерные температура, координата, время; Bi1, Bi2 – числа Био.
В безразмерном виде задача будет
; (Fo > 0; 0 < ξ < 1); (5)
; (6)
; (7)
, (8)
где
.
Будем искать решение задачи (5)–(8) в виде
, (9)
где Tc(ξ) – решение стационарной задачи с неоднородными условиями (7), (8), Θ(ξ,Fo) – решение нестационарной задачи с соответствующими однородными условиями.
Математическая постановка стационарной задачи
; (0 < ξ < 1); (10)
; (11)
. (12)
Решение задачи (10)–(12) представим в виде
, (13)
где F1, F2 – постоянные, определяемые из условий (11), (12) и имеющие вид (при T1 = T0)
,
. (14)
Аналитическое решение стационарной задачи (10)–(12), удовлетворяющее уравнению (10) и условиям (11), (12), принимает вид
. (15)
Математическая постановка нестационарной задачи
; (Fo > 0; 0 < ξ < 1); (16)
; (17)
; (18)
. (19)
Решение задачи (16)–(19) примем в виде
, (20)
где функции ϕ(Fo) и ψ(ξ) находятся из решения следующих обыкновенных дифференциальных уравнений, получаемых после подстановки (20) в (16)
; (21)
, (22)
где μ – постоянная.
Интегрируя уравнение (21), находим
, (23)
где A – постоянная интегрирования.
Граничные условия для функции ψ(ξ) находим путем подстановки (20) в (18), (19)
; (24)
. (25)
Соотношения (22), (24), (25) представляют краевую задачу Штурма – Лиувилля. Ее решение принимается в виде
, (26)
где bk(μ) – неизвестные коэффициенты.
Для определения постоянных bk(μ),
используются основные (24), (25) и ДГУ, определяемые путем выполнения искомым решением (26) уравнения (22) и соотношений, полученных после определения производных от него различного порядка, в точке ξ = 0. Любое число таких условий можно представить в виде
;
;
;
. (27)
Подставляя (26) (например, при n = 48), в граничное условие (24) и в 48 дополнительных граничных условий вида (27), для нахождения bk(μ), (k = 0,1,2,…,48) получаем систему 49 алгебраических линейных уравнений вида
(28)
Анализ системы уравнений (28) позволяет заключить, что она имеет цепочный вид с сильно разреженной матрицей коэффициентов ленточного типа. Цепочность системы состоит в том, что неизвестные b0, b1, b2, b3, …, b48 находятся из последовательного решения одного алгебраического уравнения с одним неизвестным. Положив в (28) b0 = 1, получаем следующие формулы для определения bk(μ), (k = 0,1,2,…,48):
b0 = 1;
;
;
;
;
;
;
;
;
. (29)
При известных bk(μ), (k = 0,1,2,…,48) собственные функции находятся из соотношения (26), принимающего вид
. (30)
Соотношение (30) удовлетворяет граничному условию (24), уравнению (22) и 46 соотношениям, полученным после нахождения производных от него, в точке ξ = 0. Однако оно не удовлетворяет граничному условию (25). Подставляя (30) в (25) (положив Bi1 = 3; Bi2 = 2), относительно μ получаем алгебраический полином 24-й степени. Из 24 корней полинома 6 корней выражаются действительными числами и 16 корней являются самосопряженными комплексными. Действительные корни полинома следующие:
;
;
;
;
;
;
. (31)
Точные значения первых шести собственных чисел имеют вид [14, с. 40–70]:
;
;
;
;
;
.
Учитывая, что уравнение Штурма – Лиувилля выполняется лишь при дискретных значениях μk (собственных значениях), комплексные корни, как не имеющие физического смысла, не используются. Отмечается высокая точность определения собственных чисел, отличающихся от точных их значений лишь в первом – четвертом знаке после запятой.
Для сравнения рассмотрим еще один метод определения собственных чисел, известный как метод неопределенных коэффициентов. Подставляя (26) в (22), (24) и группируя слагаемые при одинаковых степенях ξ, получаем алгебраические уравнения относительно bk(μ), (k = 0,1,2,…,48). Определяемые из их решения коэффициенты bk полностью совпадают с соотношениями (29). Собственные числа, как и выше, находятся из условия (25), и они полностью совпадают с их значениями вида (31). Следовательно, два рассмотренных метода определения собственных чисел эквивалентны. Следует, однако, отметить, что первый метод является более простым и универсальным ввиду того, что уравнение краевой задачи Штурма – Лиувилля и соотношения для производных от него различного порядка выполняются лишь в одной граничной точке ξ = 0. В связи с этим он может быть применен к решению задач с переменными физическими свойствами, с источниками теплоты и др., для которых применение второго метода затруднительно.
Докажем следующую теорему: если неизвестные функции bk соотношения (26) находить из выполнения уравнения (22) и соотношений, полученных после взятия производных от него, в граничной точке ξ = 0, то при большом числе приближений это уравнение будет выполняться во всей области изменения переменной 0 ≤ ξ ≤ 1.
Для доказательства теоремы потребуем, чтобы решение (26) удовлетворяло уравнению (22) и соотношениям, полученным после взятия k-х производных от него, в точке ξ = 0, то есть
, (k = 0,1,2,3,…). (32)
При k = 0 соотношение (33) приводится к уравнению (22). При k = 1,2,3,… оно представляет соотношения для производных различного порядка от уравнения (22). Подставляя (26) в (32), положив ξ = 0, для каждого k имеем следующие алгебраические уравнения для коэффициентов bk,
.
;
;
;
; …. (33)
Умножая (33) последовательно на ξk,
и суммируя полученные соотношения, находим
. (34)
Подставляя (26) в (22) (без привязки к конкретному значению ξ), получаем соотношение, полностью совпадающее с (34) и, следовательно, сформулированная выше теорема доказана. Отметим, что, в отличие от работы [3], здесь не требуется выполнения интеграла теплового баланса – осредненного по координате ξ уравнения (22).
Подставляя (23), (30) в (20), получаем
. (35)
Для определения постоянных An составим невязку начального условия (17) и потребуем выполнения ее ортогональности к каждой собственной функции ψi(ξ), 
. (36)
Соотношение (36) представляет систему восьми алгебраических линейных уравнений, относительно такого же количества неизвестных постоянных An.
Подставляя (15), (35) в (9), получаем решение задачи (5)–(8)
. (37)
Значения безразмерных температур, найденных по формуле (37), в сравнении с точным решением [14, с. 40–70] даны на рисунке. Видно, что для Fo ≥ 0,0001 они практически одинаковы.
Таким образом, в статье применен метод разделения переменных (метод Фурье). Получаемая при этом краевая задача Штурма – Лиувилля решается методом, предложенным авторами статьи. В ходе его реализации возникает необходимость решения систем алгебраических линейных уравнений и степенного алгебраического полинома. Их решение выполняется с использованием стандартных программных пакетов системы компьютерной алгебры MathCAD 15.0. Отсутствие в работе экспериментальных данных объясняется тем, что достоверность полученных в статье результатов подтверждена сравнением с известным точным аналитическим решением [14, с. 40–70].
Выполним исследование полученного решения для реальной конструкции керамической стенки реактора пиролиза метана при следующих исходных данных:
t0 = 20°C; t1 = 20°C; t2 = 1000°C;
δ = 0,05 м; λ = 1Вт/(м∙К); а = 1,5∙10–6 м2/с;
α1 = 60Вт/(м2∙К); α2 = 40Вт/(м2∙К).

Распределение температуры: ––––– – точное решение [14]; ○ – по формуле (37); Bi1 = 3; Bi2 = 2 Примечание: составлен авторами по результатам данного исследования
Будем рассматривать процесс теплообмена до момента времени, когда метан остается при начальной температуре, то есть прогревается только стенка.
Из анализа графиков рисунка следует, что максимальный перепад температуры на прогретой части стенки в диапазоне 0,5 ≤ ξ ≤ 1 при Fo = 0,03 (τ = 50 сек) составляет ΔΘ = 0,28 (Δt = 294°C).
Оценка температурных напряжений растяжения, найденных по приближенной формуле для тонкой стенки [15, с. 23–54]
, (38)
приводит к величине σ = 78,4 кг/мм2, что значительно превышает предел прочности для данного материала. В формуле (38) обозначено: β = 0,00001/К – коэффициент линейного расширения; Е = 20000 кг/мм2 – модуль упругости; v = 0,25 – коэффициент Пуассона.
Кроме того, на участке времени 0 ≤ Fo ≤ 0,01 отмечается высокий темп прогрева, когда температура стенки при ξ = 1 возрастает от 20°C до 176,4°C за время Δt = 37,5 сек. Столь быстрый прогрев приводит к большим значениям динамических температурных напряжений, которые также могут быть причиной появления дефектов в керамической стенке реактора пиролиза метана.
Результаты исследования и их обсуждение
Анализ результатов приводит к заключению, что ДГУ являются эффективным средством получения аналитического решения задачи теплопроводности для пластины с несимметричными условиями третьего рода. С целью теоретического обоснования полученных расчетным путем результатов приведено доказательство теоремы, по которой граничное выполнение уравнения эквивалентно его выполнению внутри рассматриваемой области. Анализируя полученные результаты, можно заключить, что научной новизной статьи является разработка методики получения высокоточного аналитического решения третьей краевой задачи, отличающееся от известных тем, что благодаря применению дополнительных граничных условий, выполняемых лишь в одной граничной точке (центре симметрии пластины), дифференциальное уравнение краевой задачи Штурма – Лиувилля выполняется во всей области определения пространственной переменной (0 ≤ ξ ≤ 1), исключая необходимость его непосредственного интегрирования.
Выводы
1. Исследования полученного решения применительно к определению температурного состояния керамической стенки реактора пиролиза метана показали, что в диапазоне начального времени температурный перепад на стенке приводит к термическим напряжениям растяжения, превышающим предел прочности для данного материала. В этом же диапазоне времени отмечается высокий темп прогрева стенки (176°C за время Δt = 37,5 сек), что может приводить к динамическим напряжениям, также превышающим предел прочности для данного материала. Для предотвращения трещинообразования в стенке реактора пиролиза метана необходимо уменьшать интенсивность подвода теплоты от греющего агента за счет постепенного повышения температуры горячих газов.
2. Данный метод может быть применен к задачам теплопроводности для цилиндра и шара, задачам с переменными физическими свойствами среды и др. Получаемые таким путем решения для цилиндрических и сферических координат не содержат функции Бесселя.