Известия РАН. Энергетика, 2019, № 4, стр. 136-148

В начальной стадии одностороннего нестационарного нагрева или охлаждения неограниченной пластины толщиной 0 ≤ х_п ≤ R_п ее температура изменяется только в поверхностном слое толщиной R, величина которого с течением времени увеличивается. В остальной части тела толщиной R_п − R температура может оставаться неизменной. Пока температурное возмущение не распространится на все тело его температурное поле принято считать полуограниченным.

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

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

Имеющиеся в литературе [1–3] аналитические решения имеют ограниченную область применения, так как не позволяют учесть реальные условия теплообмена: динамику температуры внешней среды (газов), радиационный перенос, зависимость теплофизических характеристик от температуры и др.

Численные методы расчета температурных полей более универсальны, но требуют весьма трудоемких вычислений. Для снижения трудоемкости расчетов весьма актуально создание и исследование численно-аналитических методов для моделирования температурных полей тел, некоторые из них описаны в работах [4–13].

В статьях [10, 11] был предложен метод расчета температурных полей в полуограниченных телах в форме неограниченной пластины при начальном условии Т_н(х, τ = 0) = = const, основанный на математических описаниях и алгоритмах расчета температурных полей тел простой формы, предложенных в работах [8, 9].

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

Алгоритм расчета по сравнению с [11] изменится, так как потребуется учитывать температурное поле в непрогретом слое пластины 0 ≤ х_п ≤ R_п – R и поток теплоты на поверхности х_п = 0.

Расчет процесс нагрева будет заключаться в следующем. Время нагрева пока R ≤ R_п разделяется на ряд расчетных интервалов времени Δτ. Распределение температур по сечению тела в конце каждого расчетного интервала времени Δτ аппроксимируется степенной функцией, параметры которой определяются по аналитическим выражениям, полученным решением линейной системы уравнений, которые описывают начальное и граничные условия температурного поля, а также балансы теплоты. Благодаря заранее полученным приближенным аналитическим решениям для расчетного интервала времени Δτ, решение дифференциального уравнения теплопроводности с частными производными сводится к решению обыкновенного дифференциального уравнения численным методом.

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

Рассмотрим математическое описание температурного поля ограждения или конструктивного элемента энергетической установки в виде неограниченной пластины толщиной R_п, которую в начале нагрева можно считать полуограниченным телом.

Примем, что начальное распределение температур пластины толщиной R_п, м, описывается линейной функцией

а градиенты температур на границах х_п = 0 и Х = 1 (х_п = R_п) зависят от удельных потоков теплоты где х_п − текущая координата; м; b₀, b₁− известные коэффициенты; q₀ и q − удельные потоки на границе х_п = 0 и в плоскости с координатой х = R (х_п = R_п); λ – коэффициент теплопроводности, Вт/(м ⋅ К); τ − текущее время (0 ≤ τ ≤ τ_к); τ_к − конечное время нагрева, с; R − текущая толщина поверхностного слоя (R ≤ R_п), в котором изменяется температура под действием − потока q(τ) на поверхностность х = R (х_п = R_п).

Примем, что распределение температур в поверхностном слое толщиной R в момент времени τ_{i+ 1} аппроксимируется функцией [8–10]

где Х = х/R, 0 ≤ х ≤ R, х − текущее значение координаты, изменяющейся от х = 0 при х_п = (R_п − R) до х = R при х_п = R_п (рис. 1); R − наибольшая глубина поверхностного слоя. В слое толщиной R_п − R температура не изменяется во времени пока глубина прогретого слоя R не станет равной R_п. Температурное поле Т_н(х_п, τ) при 0 ≤ τ ≤ τ_к в диапазоне 0 ≤ х_п ≤ (R_п − R) будет описываться (1)

Температура Т_н(х_п = 0, τ) и градиент температур на поверхности х_п = 0 останутся неизменными. На рис. 1 приведены обозначения координат и математическое описание распределение температур в пластине толщиной R_п для двух моментов времени.

Система с двумя текущими координатами х_п и х, где 0 ≤ х ≤ R и 0 ≤ х_п ≤ R_п, отличается от той, которая обычно применяется при описании температурного поля полуограниченного тела. Для полуограниченного тела обычно используется одна координата 0 ≤ х < ∞, где х = 0 на поверхности теплообмена.

Принятая система координат пластины позволит, не изменяя расчетные формулы использовать их для описания температурного поля пластины, когда тело прогреется на всю толщину пластины (до х_п = х = 0) и учесть теплообмен с другой стороны пластины по методу [8, 9].

Величина R будет изменяться во времени. Для конца каждого расчетного интервала значение R можно рассчитать по числу Фурье ΔFo = aΔτ_i/R²

Число Фурье следует принимать постоянным для всех расчетных моментов времени пока R < R_п, где a − температуропроводность, м/с², τ − текущее время, с. Для конца интервала i + 1 τ = τ_{i+ 1} = τ_i + ∆τ_i.

Величина числа ΔFo будет зависеть от динамики потока теплоты на поверхность тела q. Значение ΔFo должно соответствовать условию ΔFo ≡ min(R), при котором T(x = = 0, τ) = Т_н(х_п = R_п − R, τ).

Величину ΔFo можно оценить по решению дифференциального уравнения теплопроводности для полуограниченного тела. При q = const ΔFo можно принять по номограмме [3] ΔFo ≈ 0.04. В инженерных расчетах принято считать, что ΔFo ≈ 0.05.

Для расчета распределения температур в конце расчетного интервала времени τ_{i+ 1} = = τ_i + ∆τ_i по формуле (4) необходимо определить коэффициенты а₀, а₁, а₂ и показатель степени n.

Формулы для вычисления коэффициентов a₀, a₁, а₂ и показателя степени n можно получить, используя следующие уравнения.

1. Граничное условие, уравнение (3), которое для конца интервала времени при х = R и Х = 1 с учетом выражения (4) примет вид

где q_к − поток теплоты в конце расчетного интервала времени q_к(τ_{i+ 1}).

Из (3а) можно выразить коэффициент n или a₂ для момента времени τ_{i+ 1}

2. Уравнение баланса теплоты с учетом (1) в плоскости Х = 0 и х_п = R_п − R_{i + 1}

Потоки теплоты в (9) с учетом (1) и (4) можно записать в виде

и получить из (9) с учетом (1) и (4) формулу для расчета а₁

3. Уравнение равенства температур в плоскости Х = 0 и х_п = R_п − R_{i + 1}

температуры для которого выразим из (1) и (4) и получим из (13), (14) и (15) формулу для расчета а₀ ${\text{\;}}\left( {14} \right)$

4. Уравнение баланса теплоты для прогретого слоя толщиной R для интервала времени Δτ_{i+ 1}в конечных разностях, которое запишем в следующем виде:

где c·− объемная теплоемкость, Дж/(м³ · К); q_ср·− среднеинтегральный поток теплоты для расчетного интервала времени; T_{ср, ни} и Т_{ср,i+ 1} – среднемассовые температуры в начале и конце интервала времени Δτ_{i+ 1} = τ_{i+ 1} − τ_i для слоя толщиной R_{i + 1}.

Величина T_{ср, ни} в начале расчетного интервала времени τ_i должна определяться по температуре для конца предыдущего интервала времени Т_ср,i по формуле

Обратим внимание, что T_ср,i − среднемассовая температура в конце i-го интервала времени относится к слою R_i, а среднемассовая температура в начале i + 1 интервала времени T_ср, ни − к расчетному слою R_{i + 1} (см. уравнение баланса теплоты (17)). Поэтому на i + 1 интервале времени значение среднемассовой температуры в начале i + 1 интервала времени T_ср, ни должно определяться с учетом T_ср, i и температуры непрогретой части слоя (R_{i + 1} − R_i), которая равна T_{н, ср,i+ 1}.

Выражение для расчета среднемассовой температуры непрогретого слоя 0 ≤ х_п ≤ R_п − R_i в конце интервала Т_{н, ср, i+ 1} можно получить интегрированием (5)

где ${{x}_{{{\text{п}},i}}} - {{x}_{{{\text{п}},\,i + 1}}} = {{R}_{{i + 1}}} - {{R}_{i}}.$

Значение среднемассовой температуры слоя R_i + 1 в конце интервала Т_{ср, i + 1} определится интегрированием (4) [8, 9]

С учетом (20) уравнение баланса теплоты (17) запишется в виде

где средний поток теплоты q_{ср, i+ 1} определится по потокам теплоты в начале q_н, i + 1 и конце q_{к,i+ 1} интервала времени

Из уравнения (21) выразим коэффициент а₂

где

С учетом выражения (7) формулу (24) для расчета n представим в виде

и запишем решение квадратного уравнения (25) относительно а₂

Отметим, что при ΔТ < 0 коэффициент а₂ следует рассчитывать по формуле (8), задавая n = 4.

Зная коэффициент а₂, показатель степени n можно определить по формуле (7).

Таким образом, по довольно простым формулам (12), (16) и (26) можно определить а₁, а₀, а₂ и n по формуле (7), рассчитать распределение температур в конце расчетного интервала времени τ_{i+ 1} по (4).

Сравнительно простые формулы (12), (16), (26) и (7) − это аналитическое решение дифференциального уравнения теплопроводности при граничных условиях второго рода для интервала времени Δτ_{i+ 1} при известных значениях T_{ср, ни} и Т_н(х_п = R_п − R).

Расчет температур поверхности Т_{1,i+ 1} = Т(Х = 1, τ_{i+ 1}) и среднемассовой температуры в конце интервала времени Т_{ср,i+ 1} можно выполнить по формулам

которые получены из (4) при Х = 1 и из выражения (20).

Для первого интервала среднемассовая температура T_{ср, ни} в слоя толщиной R₁ начале интервала (i = 0) определится из начальных условий по формуле (19) при R_{i = 0} = 0 и x_п,i = R_п − R₀ = R_п.

Для последующих интервалов i > 1 значение T_{ср, ни} найдется с учетом изменения размера прогретого слоя по формуле (18) с учетом (19).

Расчет распределений температур в последующие моменты времени τ_{i+ 1} = τ_i + Δτ можно выполнять с помощью известных процедур решения обыкновенных дифференциальных уравнений, например, методов Эйлера, Эйлера-Коши, Рунге-Кутта, ломаных и др.

Для заданной функции q(τ), расчет по описанным выше формулам не вызовет серьезных затруднений, так как будет проводиться по простым алгебраическим выражениям. Однако для сопряженной задачи теплообмена система уравнений для расчета а₁, а₀, а₂ и n будет нелинейной.

Обратим внимание, что в формулу (26) для расчета а₂ с учетом (24) входят q₀(х_п) и средний поток теплоты q_{ср,i+ 1}, который можно определить по потокам теплоты в начале q_{н,i+ 1} и конце q_{к,i+ 1} интервала времени по (22).

Величина q₀(х_п) зависит от начальных условий и легко находится по (10).

В начале расчета температур на интервале Δτ_{i+ 1} имеются известные данные только для вычисления потока теплоты в начале интервала q_{н,i+ 1} по известной температуре поверхности Т_1,i. Температура поверхности Х = 1 Т_{1,i+ 1} = Т(Х = 1, τ_{i+ 1}) и величина потока теплоты в конце интервала времени q_{к,i+ 1} неизвестны, поэтому их значения приходиться принимать приближенно. При расчете по методу Эйлера принимается, например, условие q_{ср,i+ 1} ≈ q_{к,i+ 1} ≈ q_{н,i+ 1}.

Для материалов с низкой теплопроводностью (кирпича, тепловой изоляции) и металлов при интенсивном теплообмене характерна высокая скорость роста температуры обогреваемой поверхности, которая приводит к резкому изменению потока теплоты в пределах расчетного интервала времени. В таких условиях задание q_{ср,i+ 1} ≈ q_{к,i+ 1} ≈ q_{н,i+ 1} вызывает значительные погрешности расчета или приводит к абсурдным результатам.

В работах [8, 10] предложена оригинальная методика для оценки значений температуры поверхности (Х = 1) Т_{1,i+ 1} и q_{к,i+ 1}. Методика основана на допущении постоянства в пределах Δτ_{i+ 1} приведенного коэффициента теплообмена.

Модифицируем ее для данного математического описания.

Рассмотрим процедуру получения первого приближения температуры поверхности $Т_{{{\text{1}},\,i + {\text{1}}}}^{{\text{п}}}.$

Поток теплоты q на поверхность Х = 1 в расчетные моменты времени в общем случае может зависеть от температур газов Т_г(τ), коэффициентов радиационного σ, Вт/(м² · К⁴), и конвективного α_к, Вт/(м² · К), теплообмена.

Запишем выражение для приведенного коэффициента теплообмена в начале расчетного интервала

где q_{н,i+ 1} − радиационно-конвективный поток теплоты, Вт/м², вычисленный по известной температуре поверхности T_{1,i+ 1}в начале i + 1 интервала, равной Т_1,i,

Выразим первое приближение для потока теплоты в конце i+1 интервала через приведенный коэффициент теплообмена и неизвестную пока температуру поверхности в конце интервала $Т_{{{\text{1}},\,i + {\text{1}}}}^{{\text{п}}}$

Используя выражения (27) и (28) запишем разность температур $Т_{{{\text{1}},\,i + {\text{1}}}}^{{\text{п}}}$ и Т_{ср,i+ 1} с учетом формулы(8) для а₂

Из уравнения баланса теплоты (17) получим формулу для Т_{ср,i+ 1} и подставим ее в (32). Преобразовав (32) с учетом (17) получим выражение для расчета первого приближения температуры $Т_{{{\text{1}},\,i + {\text{1}}}}^{{\text{п}}},$ которую используем для расчета первого приближение для потока теплоты в конце i + 1 интервала по формуле (31):

где $w = \frac{{\Delta {{{\tau }}_{{i + 1}}}}}{{2c{{R}_{{i + 1}}}}}.$ Единица в показателе степени обозначает первое приближение.

Величина $Т_{{{\text{1}},\,i + {\text{1}}}}^{{\text{п}}}$ является приближенной, так как коэффициент α_{пр, н} справедлив не для всего интервала, а только для его начала. Значение показателя степени n^п для первого интервала времени следует принять, например n^п = 3. Для второго и последующих этапов его величину можно принять по предыдущему интервалу времени. Исследования показали, что задание n^п для первого интервала времени в диапазоне 3 < n^п < 6 и 0.03 < Fo < 0.07 практически не отражается на результатах расчета температуры поверхности тела Т(1, τ).

По известной температуре $Т_{{{\text{1}},\,i + {\text{1}}}}^{{\text{п}}}$ можно рассчитать приближенное значение потока теплоты в конце интервала q_{к,i+ 1} по формуле (31), используя $Т_{{{\text{1}},\,i + {\text{1}}}}^{{\text{п}}}$ и другие параметры теплообмена соответствующие τ_{i+ 1}. По q_{н,i+ 1} и q_{к,i+ 1} можно вычислить средний для интервала времени поток теплоты q_{ср,i+ 1} (см. (22)). Далее определить а₀, а₁ и показатель степени n и по функции (4), рассчитать распределение температур в конце расчетного интервала времени T(X) или Т_{1,i+ 1}. и Т_{ср,i+ 1} по (27), (28).

В некоторых случаях значение Т_{1,i+ 1} можно уточнять простыми итерациями, используя его в качестве следующего приближения вместо T_{1, н,i+ 1} при расчете α_пр по (20).

Для снижения трудоемкости расчета от расчета $Т_{{{\text{1}},\,i + {\text{1}}}}^{{\text{п}}}$ можно отказаться, когда потоки теплоты незначительно различаются в пределах Δτ, то есть когда q_{ср,i+ 1} ≈ q_{н,i+ 1} ≈ ≈ q_{к,i+ 1}. Это возможно для металлов или для тел с низкой теплопроводностью, когда слабо изменяется интенсивность теплообмена.

Рассмотрим пример расчета нагрева материала пластины (с = 1.5 × 10⁶; λ = 0.8; a = = 5.333 × 10^–7; R = 0.2) со стороны x_п = R_п = 0.2 в среде Т_г = 600, α = 60 при начальном распределении температур Т_н(x_п) = 400–500x_п (b₀ = 400, b₁ = –500) (функция Т_н(x_п) показана на рис. 1).

Определим параметры температурного поля в конце i = 1 периода нагрева.

Примем i = 0; n¹ = 3; R_{i = 0} = 0; Fо = 0.05; Δτ₁ = τ₁ = 10; Т_г = 600.

Рассчитаем глубину прогрева по формуле (6)

и коэффициент a₁ = b₁R₁ = −500 × 0.010332 = −5.164.

Определим диапазон интегрирования для расчета средней температуры Т_{ср, ни} слоя R₁ = 0.01033 в начале интервала времени x_{п, 0} = 0.2; x_{п, 1} = 0.2 − 0.01033 = 0.1897 и ее значение

Потоки теплоты в x_{п, 1}q₀ = λb₁ = 0.8(–500) = –400 и на поверхности x = R (X = 1) в начале интервала q_н = 60(600 – 300) = 18000.

Приведенный коэффициент теплообмена α_пр = α = 60 и комплекс w₁

Первое приближение для температуры поверхности в конце интервала времени

Первое приближение потока теплоты в конце первого интервала Δτ

Средний поток теплоты на интервале Δτ₁

Коэффициент а₀ = 400 + (–500) × 0.1897 = 305.16.

Величина ΔТ по (24)

Коэффициент а₂ по (26)

Показатель степени n по формуле (7)

Температура поверхности Т_{1,i+ 1} (x = R, X = 1 и x_п = R_п)

(Значение Т_{1,i+ 1} отличается от первого приближения на 358.7 – 351.1 = 7.6 К.)

Средние температуры прогретого слоя ${{R}_{1}}$ в конце интервала Δτ

и непрогретого слоя и всей пластины

Параметры в начале второго интервала τ₂ = 20, Δτ₂ = 20 – 10 = 10, i = 1

Описанная процедура повторяется для следующих интервалов времени, пока R < R_п.

Результаты расчета для времени 0 ≤ τ ≤ 3600 с, пока R < R_п, приведены на рис. 2–4.

На рис. 2 показана динамика температур газа, поверхности х_п = R_п, средних температур прогретого слоя и всей пластины.

На рис. 3. Показана динамика изменения глубины R прогрева пластины.

По рис. 2 и 3 можно сделать вывод, что для принятых условий, после прогрева пластины наполовину толщины R > R_п и τ > 900 с средние температуры прогретого слоя и всей пластины становятся практически равными.

На рис. 4 показаны распределения температур по толщине пластины, рассчитанные по формулам (1) и (5) для нескольких моментов времени.

Рассмотрим результаты расчета радиационно-конвективного нагрева (α = 40, σ = = 4 × 10^–8, Т_г(τ) = 900 = const) стальной пластины R_п = 0.2, свойства которой зависят от температуры (коэффициент теплопроводности λ = 63.41 – 0.03256Т, температуропроводность а = 18.1 × 10^–6 – 1.34 × 10^–8Т). Начальные распределения температур Т_н(x_п) = b₀ + b₁x_п задавались для двух вариантов:

• 1 – b₀ = 400, b₁ = –500 (внешний нагрев со стороны х_п = 0);

• 2 – b₀ = 300, b₁ = +500 (внешнее охлаждение со стороны х_п = 0).

Величина ΔFo принималась равной 0.051, а Δτ = 15 с.

Обратим внимание, что величины а в (6) и λ в (3), (7), (8) рассчитывались по Т(Х = 1, τ), а λ в (10) – по температуре Т_н(x_п = R_п – R).

На рис. 5 показана динамика температур поверхности х_п = R_п, средней температуры всей пластины для двух вариантов задания начальных условий.

На рис. 5 видно, начальное распределение температур оказывает значительное влияние на скорость роста температуры поверхности.

При внешнем нагреве со стороны х_п = 0 Т₁ увеличилась за 150 с на 70 К, а при охлаждении – только на 20 К.

При тех же условиях, но при температуре газа Т_г(τ) = 600 = const, температуры Т₁ и Т_ср пластины даже снижались при нагреве с верху (х_п = R_п) и охлаждении с низу (х_п = 0). В этом случае величина ΔТ < 0 (см. (25)) и расчет а₂ выполнялся по формуле (8) и n = 4. При n = 3 и n = 5 погрешность расчета увеличивалась, но не более чем на 5 К.

Отметим, что при тех же условиях теплообмена были выполнены расчеты методом конечных разностей. Расхождений результатов не превышало 1%.