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

Однако следует отметить, что, применительно именно к оценке теплового режима зданий, до настоящего времени основное внимание уделялось так называемому «регулярному режиму», когда тепловым возмущением охвачена вся толща материала конструкции, особенно при перерывах в теплоподаче [1], либо периодическим колебаниям [2] (теория теплоустойчивости).

Однако при автоматизации климатических систем основное значение имеет именно скорость распространения температурных волн в массиве ограждений, и здесь требуется несколько иной математический аппарат. В этой области также существует ряд работ, например [3, 4], где задача рассматривается, тем не менее, с такими же упрощающими допущениями, как для регулярного режима. В последние годы возникают также некоторые публикации, в которых особое внимание обращается на более тщательный учёт граничных условий при численном моделировании исследуемых процессов [5, 6], либо содержащие результаты экспериментальных исследований переходных явлений в ограждающих конструкциях [7], а также посвящённые использованию экономических методов для комплексного изучения условий функционирования оболочки здания [8]. Однако достигаемые при этом результаты сложно назвать доступными для практики проектирования.

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

здесь x — текущая пространственная координата, соответствующая расстоянию от наружной поверхности стенки до исследуемого сечения, м; t — температура материала стенки в этом сечении в момент времени τ [с] с начала разогрева или охлаждения, К; a = λ/(сρ) — коэффициент температуропроводности материала стенки, м²/с; λ — его теплопроводность, Вт/(м·К); с и ρ — удельная массовая теплоёмкость [Дж/(кг·К)] и плотность [кг/м³], соответственно.

Чтобы обозначить необходимое нам постоянство плотности теплового потока через поверхность q [Вт/м²] при x = 0, нужно записать граничное условие второго рода:

Второе условие, очевидно, формулируется как t = 0 при x → ∞, если для простоты рассматривать избыточную температуру по отношению к начальной.

Существующие решения рассматриваемой задачи, в частности [9], предусматривающие запись уравнения в форме (1) относительно величины q, что позволяет добиться автомодельности по отношению к безразмерному комплексу x/(2КОРЕНЬ[aτ]), и последующее вторичное интегрирование после подстановки результата в (2), показывают, что интересующая нас в наибольшей степени температура t0 на поверхности стенки (при x = 0) пропорциональна произведению КОРЕНЬ[τ] на комплекс q/КОРЕНЬ[λcρ].

К данному выводу можно прийти и при использовании метода анализа размерностей, но аналитическое решение даёт и конкретную величину числового коэффициента, который в данном случае оказывается равным 2/КОРЕНЬ[π] ≈ 1,13.

Однако некоторые экспериментальные данные и результаты численных расчётов, полученные при участии автора [10, 11], свидетельствуют в пользу иного значения коэффициента, приближающегося к КОРЕНЬ[π] ≈ 1,77.

Для разрешения данного противоречия можно воспользоваться оценочным подходом, основанным на получении приближённого решения методом Ритца [12]. Варьируем форму профиля температуры в зоне распространения температурной волны при точном выполнении граничных условий на поверхности стенки и на фронте теплового возмущения, считая коэффициент пропорциональности k для глубины его проникновения ∆ = kКОРЕНЬ[aτ] равным 2КОРЕНЬ[π] ≈ 3,55 в соответствии с фундаментальным решением дифференциального уравнения теплопроводности [9], а также с результатами численных расчётов [11]. В качестве изменяемого параметра будет выступать показатель степени n в полиномиальной аппроксимации температурного профиля. Тогда, если записать такое приближенное выражение для температуры в общем виде, получим:

t = t0(1 — z)n, (3)

где z = x/∆ = x/(kКОРЕНЬ[aτ]) — относительная безразмерная координата.

Подставляя данную зависимость в условие (2), находим соотношение для t0:

то есть поверхностная температура в любом случае будет пропорциональна КОРЕНЬ[τ], а искомый коэффициент пропорциональности для принятой формы записи решения оказывается равным k/n.

Если теперь учесть, что все остальные параметры, входящие в формулу для t0, постоянны и не связаны ни с x, ни с τ, выражение для невязки по методу Ритца можно получить в упрощённом виде, для условной относительной температуры

Поскольку КОРЕНЬ[τ] не зависит от n, находим минимум интеграла от квадрата невязки φ(z, n) в естественных пределах z от 0 до 1. Из-за сложности полученного выражения это проще всего сделать численными методами, и тогда оказывается, что оптимальное значение n = 2,88. В то же время общее количество теплоты [Дж/м²], поглощённое единицей площади стенки за время τ, очевидно, будет равно:

где m — коэффициент, показывающий, во сколько раз средняя величина температуры в зоне распространения температурной волны меньше, чем t0. Нетрудно получить, что при полиномиальной аппроксимации типа (3) m = n + 1. Теперь для произвольного k с учётом (4) и выражения для ∆ находим:

Однако по физическому смыслу данное произведение должно равняться просто qτ. Следовательно, для наилучшего приближения имеем k2/[n(n + 1)] = 1, откуда при k = 2КОРЕНЬ[τ] вычисляем:

n = 0,5×(КОРЕНЬ[16π + 1] — 1) ≈ 3,08,

что достаточно близко к данным предыдущего расчёта.

На рис. 1 приведён соответствующий график безразмерной относительной температуры t/t0: красной линией изображены результаты численного решения, а голубой линией — профиль температуры при n = 3,08. Синяя линия соответствует профилю температуры при n = 2,88.


Рис. 1. Безразмерные профили температуры в пределах зоны распространения температурной волны (результаты численного расчёта и аппроксимаций)

Численное решение выполнено с использованием составленной автором программы для ЭВМ на языке Fortran, основанной на конечно-разностной аппроксимации уравнения (1) по чисто неявной схеме [12, 13], для нескольких вариантов теплотехнических свойств материала и толщины стенки, и все полученные графики в пределах толщины линии накладываются друг на друга, поэтому можно считать его достоверным. Отсюда же подтверждаем равенство k = 2КОРЕНЬ[τ], поскольку видно, что вблизи поверхности стенки при n = 3,08 совпадение с приближением (3) при использовании данного значения k очень хорошее, то есть граничное условие здесь действительно выполняется точно.

В то же время для x = ∆ величина t, даваемая приближенным выражением, строго равняется нулю, то есть соблюдается и условие на фронте температурной волны, но это следует из самой формы записи (3). При этом числовой коэффициент k/n в формуле для t0 оказывается равным 2КОРЕНЬ[τ]/3,08 ≈ 1,15, что практически не отличается от величины 2/КОРЕНЬ[τ] из решения [9]. В случае n = 2,88 мы имеем наименьшее отклонение в среднем по глубине проникновения температурной волны — как, собственно, и должно быть при минимизации невязки в целом по отрезку, а k/n = 1,23, однако вблизи x = 0 расхождение с численным решением несколько больше.

Тем не менее, при осуществлении аппроксимации (1) по явной схеме с максимально возможным из условия устойчивости шагом по времени, равным h2/(2a), где h — шаг по пространственной координате [м], что позволяет предельно упростить и ускорить расчёт, результаты получаются несколько иными. В этом случае на каждом следующем, (j + 1)-м временнóм шаге температура в i-м узле сетки рассчитывается как полусумма значений в обоих соседних узлах в предыдущий, j-й момент времени [12, 13]:

Тогда данные вычислений, выполненных опять-таки при различных сочетаниях характеристик материалов и размеров конструкции, свидетельствуют о том, что при измельчении шага значение k/n стремится к 1,592, что уже приближается к полученному в [10, 11]. На рис. 2 рассчитанный таким способом безразмерный профиль t/t0 показан красной линией, голубым цветом отмечена его аппроксимация по типу (3) для n = 2,15, что отвечает максимальному приближению граничного условия на поверхности, а синяя линия изображает профиль (3) при n = 1,9, что означает аппроксимацию в среднем.


Рис. 2. Безразмерные профили температуры в пределах зоны распространения температурной волны (результаты численного расчёта и аппроксимаций)

Таким образом, вид распределения температуры в прогретой (охлаждённой) зоне и, соответственно, коэффициенты пропорциональности при t0 оказываются существенно зависящими от применяемого метода расчёта, чего, вообще говоря, быть не должно. Поэтому для окончательного разрешения этого вопроса было бы целесообразно выполнить некоторые дополнительные исследования теоретического и экспериментального характера, в том числе численные и натурные. В этом случае дальнейшие результаты, описанные автором в работе [11], в которой рассматриваемый коэффициент пропорциональности для упрощения записи и введения некоторого запаса на несовпадение величины t0 и температуры внутреннего воздуха в помещении tв был принят равным 2,0, получат более полное подтверждение и возможность применения выявленных зависимостей в инженерной практике.