CC BY
25
УДК 536.2; 519.633.2
Восстановление коэффициентов теплопроводности и объемной теплоемкости как функций температуры по экспериментальным данным установки для исследования рыхлых грунтовых пород
Д-р техн. наук Е. С. ПЛАТУНОВ, С. В. СЕРГЕЕВ
Санкт-Петербургский государственный университет низкотемпературных и пищевых технологий
191002, Санкт-Петербург, ул. Ломоносова, 9
The existing analytical method to treat experiments for defining thermal physical characteristics of moisture containing materials as functions of temperature on the basis of a dynamic method of free cooling results in incorrect heat conductivity coefficients in the range of phase transfer temperatures. A numerical method to estimate reliability of coefficients computed by a proposed method is presented. To define the values of recovered thermal physical characteristics more precisely a method to treat experimental data is proposed; it is based on a numerical solution of the coefficient inverse problem of heat conductivity where data logged in the experiment proper are used as boundary condition, namely temperature inside the specimen and heat flow on its surface.
Key words: monotonous regime, moisture containing materials, inverse problem, heat measuring cell, dynamic method, heat conductivity, heat capacity.
Ключевые слова: монотонный режим, влагосодержащие материалы, обратная задача, теплоизмерительная ячейка, динамический метод, теплопроводность, теплоемкость.
Для определения теплофизических характеристик (ТФХ) влагосодержащих материалов используют динамические методы свободного охлаждения. Экспериментальная часть метода представляет собой следующее: теплоизмерительная ячейка с образцом цилиндрической формы внутри помещается в холодильную камеру, где с интервалом в 5—10 с регистрируется температура в центре и на его поверхности. В результате получаем два массива экспериментальных данных температур центра и поверхности. После сглаживания имеем две экспериментальные кривые «температура—время» (рис. 1).
Кроме того, в установке, реализующей данный метод, есть возможность измерять тепловой поток, проходящий через поверхность образца. Для него также получаем массив данных, который при сглаживании дает экспериментальную кривую «тепловой поток—время» (рис. 2).
Режим свободного охлаждения образца особенно подходит для исследования влагосодержащих материалов, поскольку оказывается очень близким к монотонному [1). Это обеспечивается как самой конструкцией установки, так и наличием фазового перехода в образце,
который, с одной стороны, помогает выдерживать требуемый по условиям монотонности перепад температур внутри образца, а с другой — является «проблемной» областью для расчета ТФХ. Математическое описание монотонного режима подразумевает целый ряд оправданных допущений, которые приводят к достаточно простым и компактным аналитическим соотношениям для непосредственного вычисления значений ТФХ.
Тем не менее, приходится вводить определенные поправки в конкретные соотношения с целью учета ряда факторов. К примеру, конечные размеры образца приводят к искажению его температурного поля за счет теплообмена торцевой поверхности с внешней средой, в то время как соотношения подразумевают одномерное температурное поле внутри образца. Относительно большие перепады температуры внутри образца, которые на определенных этапах эксперимента не удается выдерживать в требуемых пределах, также требуют внесения поправок. Величины указанных поправок могут быть оценены специальным образом.
Т,°С
Рис. 1. Характерный вид экспериментальных температурных кривых:
1 — температура центра; 2 — температура поверхности
Вычисление значений ТФХ производится для каждого момента времени из массива экспериментальных данных. В результате получается расчетный массив дискретных значений ТФХ, каждому из которых ставится в соответствие температура образца, которая наблюдалась в соответствующий момент времени. При сглаживании данного массива получаются кривые ТФХ как функции температуры. Следует отметить, что в данном случае имеется в виду среднеобъемная температура образца, которая представляет собой среднеарифметическое между температурой в центре образца и на поверхности.
Анализ кривых ТФХ как функций температуры, рассчитанных указанным методом, выявляет физически необъяснимое поведение, в частности, кривой коэффициента теплопроводности в диапазоне температур фазового перехода. Причем подобная ситуация наблюдалась при исследовании различных образцов. Это вызвано тем, что затвердевание образца происходит в некотором интервале температур, где зависимость его искомых ТФХ от температуры становится сугубо нелинейной, а соотношения монотонного режима утрачивают силу.
Для оценки достоверности рассчитанных выше ТФХ предлагается следующий математический эксперимент. Смоделируем работу тепловой ячейки установки в виде двумерной нелинейной задачи теплопроводности с известными из эксперимента граничными условиями (тепловой поток как функция времени). При этом ТФХ как функции температуры нам также уже известны. Методом конечных разностей на подробной пространственно-временной сетке рассчитываем поле температур внутри образца как функцию времени. Далее выделяем изтемпе-ратурного поля две точки (в центре образца и на границе) и получим две кривые температура—время. Затем сравни-
£>, Вт
Рис. 2. Характерный вид кривой теплового потока на поверхности образца
ваем полученные расчетным путем температурные кривые с экспериментальными температурными кривыми, приведенными на рис. 1. Если используемые при решении прямой задачи зависимости А (Г) и с (Т) соответствуют истине, то температурные кривые должны совпадать в рамках погрешности численного решения и погрешности измерения температур в опыте.
Такой математический эксперимент был проведен с тремя различными образцами грунта. Во всех случаях наблюдалось значительное расхождение рассчитанных температурных кривых с экспериментальными данными (рис. 3).
Г, 0 С
Рис. 3. Сопоставление расчетных и экспериментальных температурных кривых
Такое расхождение кривых свидетельствует о том, что аналитические соотношения для расчета коэффициентов теплоемкости и теплопроводности приводят к некорректным результатам именно в диапазоне температур фазового перехода, в то время как в диапазоне температур однофазной области они дают достоверный результат.
В данной статье предлагается другой способ обработки экспериментальных данных, приводящий к получению более точных и физически корректных результатов в зоне температур фазового перехода, основанный на численном решении обратной задачи теплопроводности. Математическая модель процесса охлаждения образца в теплоизмерительной ячейке представляется двумерным нелинейным уравнением теплопроводности с известными из эксперимента граничными условиями: тепловыми потоками и температурами на поверхности и в центре образца. При этом тепловой поток в центре образца равен нулю в силу осесимметричности задачи, атепловой поток на поверхности и обе температуры являются функциями времени. Таким образом, имеем коэффициентную обратную задачу, когда по известным температурам в двух точках образца и тепловому потоку на поверхности находятся коэффициенты теплопроводности и объемной теплоемкости как функции температуры [2]. Теория и методы решения обратных задач развиты достаточно хорошо, и это позволяет применять их для построения алгоритма расчета ТФХ, используя в качестве исходных данных кривые температур, показанные на рис. 1.
Учтем, что в теплоизмерительной ячейке [3] образец влагосодержащего материала имеет форму ограниченного цилиндра с высотой Я, радиусом Д и начальной температурой Т0. В начале опыта (т > 0) образец начинает охлаждаться, через его основание и боковую поверхность проходит тепловой поток плотностью д = /(т), при этом температура боковой поверхности (стенки) Тст и в центре образцаТц изменяется стечением времени. Функции Тст(т) и Тц(т) измеряются в процессе эксперимента. Задача состоит в том, чтобы с их помощью установить ТФХ исследуемого образца: теплопроводность А(Т) и объемную теплоемкость с(Т), являющиеся функциями температуры.
Сформулированная задача сводится к решению квазилинейного двумерного уравнения теплопроводности [4]:
,тдТ 1 д (х/т (дТ дТ\\
С( ( ^ (ж + ж)) +
д Л/^ч (дт + — ( ЧТ)
дТ
дії
<% V 4 ' \дг с начальным условием
Т(г,/і,0) = Т0
и граничными условиями
<7о/0
й^дТ сГг ~дт
(1)
(2)
дТ
дг
— О,
г=0
дТ
Ж
дТ
'Ж
/г=0
дТ
= ^’ Ж
= 0;
(3)
Н=Н
= Чт(т),
~=Я
Т(0,Я/2,т) = Т*(т), Г(Я,Я/2,т) = Тс*т(т). (4)
Здесь Т(г, Л, 0) — температура в точке с координатами г, И в момент времени т;
д0 — скрытая теплота фазового перехода;
6(Т) — безразмерный параметр состояния (при в = 0 образец полностью разморожен, при в = 1 образец полностью заморожен);
Чг{т) = <27і(т) — плотность теплового потока на боковой и торцевой поверхностях образца.
Величина теплового потока может быть измерена на градиентном тепломере в каждый интересующий нас момент времени. При этом тепловой поток представляет собой плавную функцию времени, сначала возрастающую, затем асимптотически убывающую.
Задача восстановления функций \(Т) и с(Т) из условий (1—4) сводится к минимизации функционала вида
Т гп
I(А, с) = I { [Т(0,Я/2,т)-Тц’(0,Я/2,т)]2 + о
+ [Т(Я,Я/2,т) -Т*(Я,Я/2,т)]2 }с1т -* шт < 52. (5)
Эта процедура требует многократного решения квазилинейного уравнения теплопроводности (1) с условиями (2—4). Конкретное решение уравнения (1) осуществляется численным методом конечных разностей.
При переходе к конечным разностям дифференциальное уравнение (1) и условия (2—4) образуют систему квазилинейных уравнений, которая решается на каждом временном шаге с помощью итерационного процесса и метода прогонки.
Минимизацию функционала (5) наиболее естественно осуществлять градиентным методом [5]. В свою очередь, определение градиента функционала квадратичной невязки (5) сводится к решению системы разностных уравнений прямой задачи (1—4), решению сопряженной системы уравнений и вычислению компонент градиента.
Таким образом, процедура поиска функций А(Т) и с(Т) сводится к нахождению узловых точек А(АЬ А2,..., А„) и с(с1, сг,..., с„), число которых определяется шагом выбранной сетки, и к последующему способу их интерполяции.
Апробирование предложенного метода обработки опыта производилось на модельной задаче. В частности, важным оказался выбор способа минимизации функционала (5). Для его минимизации первоначально использовались три метода: градиентный метод, метод наискорейшего спуска и метод сопряженных градиентов. Лучшие
результаты дал метод сопряженных градиентов. При этом анализ показал, что во всех трех случаях все же не удается обеспечить равномерную сходимость решения у границ исследуемой области. Причем улучшение сходимости не наблюдается даже при использовании интерполяции между узловыми точками более высокого порядка, чем линейная, в частности квадратичная, кубическая, сплайн-интерполяция. Такая ситуация имела место при восстановлении ТФХ поданным реальных опытов с глиной, яблоками и мясом, проведенных Тамбулатовой Е. В. (см. рис. 4, 5). Это может быть связано с тем, что компоненты градиента функционала (5), соответствующие граничным узлам функции А(АЬ А2,..., А„), в несколько раз меньше компонент, соответствующих внутренним узлам. Последнее обусловлено характером температурного поля в образце, системой исходных данных и величиной шага табличного представления функций А(Т) и с(Т). В результате при минимизации функционала граничные точки области движутся к искомому решению значительно медленнее остальных, что и вызывает искажение решения у границ.
Наиболее эффективным устранением неравномерной сходимости оказалось использование процесса последовательной минимизации, состоящего в последовательном увеличении числа узлов интерполяции. Первона-
X., Вт/(м К)
Рис. 4. Коэффициент теплопроводности яблока: Я — аналитическая обработка данных;
□ — решение обратной задачи при числе узлов интерполяции N = 8;
Д — решение обратной задачи при числе узлов интерполяции N = 16
чальная минимизация на крупной сетке обеспечивает хорошее начальное приближение для последующей минимизации на более мелкой сетке. При таком подходе заметно улучшается сходимость решения у границ области, а также в «местах» быстрого изменения искомой функции, которое наблюдается в диапазоне температур фазового перехода.
Следует также отметить, что простое измельчение сетки хотя и приводит к улучшению сходимости, но значительно увеличивает время счета и требует больших ресурсов оперативной памяти компьютера. Это вызвано многократным решением задачи по нахождению температурного поля. При этом в памяти должен храниться тензор рассчитанных температур высокого ранга.
Таким образом, приведенный способ обработки экспериментальных данных имеет явное преимущество при использовании его в области сильной зависимости коэффициентов переноса от температуры, а в окрестностях границ исследуемой области температур этот способ дает менее точные результаты, чем аналитический. На рис. 4 видно, что даже при расчете на самой грубой сетке рассчитанные значения коэффициента теплопроводности образуют плавную кривую на всей исследуемой области температур.
С, кДж/(м3-К)
Рис. 5. Коэффициент объемной теплоемкости яблока: • — приближенный аналитический расчет;
□ — результат расчета при N = 7;
Д — результат при N — 15
Описанный численный алгоритм восстановления коэффициентов А (Г) и с(Т) оказывается в целом более точным в сравнении с аналитическим методом расчета при условии использования алгоритма последовательной минимизации. Однако существенное ухудшение сходимости численного алгоритма у границ исследуемой температурной области ставит под вопрос его универсальность. В этом мы убедились, используя найденные в обоих случаях коэффициенты для решения прямой задачи теплопроводности на цилиндрических образцах, используемых в теплоизмерительной ячейке. Полученные временные зависимости температур в соответствующих точках модели сравнивались с экспериментальными данными. Использование коэффициентов, восстановленных численно, дало достаточно хорошее совпадение кривых температур на всей временной области, кроме окрестности границ. При использовании же коэффициентов, полученных аналитически, наблюдалось значительное расхождение кривых температур в области сильной зависимости ТФХ от температуры.
Список литературы
1. Платунов Е. С., Баранов И. В., Буравой С. Е., Курении В. В. Теплофизические измерения / Под общ. ред. Е. С. Платунова. - СПб.: СПбГУНиПТ, 2010.
2. Денисов А. М. Введение в теорию обратных задач: Учеб. пособие. — М.: МГУ, 1994.
3. ЕІлатунов Е. С., Баранов И. В., Куслиева Е. В. Автоматизированный прибор для измерений теплофизических характеристик влагосодержащих материалов // Вестник МАХ. 2009. № 3.
4. Актуальные проблемы механики, прочности и теплопроводности при низких температурах. Моделирование технологий замораживания грунтов искусственным холодом: Материалы IX научно-технической конференции / Под ред. В. В. Улитина. — СПб.: СПбГУНиПТ, 2003.
5. Карчевский А. Л. Корректная схема действий при численном решении обратной задачи оптимизационным методом // Сибирские электронные математические известия. 2008. Т. 5.
