Введение в моделирование нестационарной теплопередачи в твёрдых телах с помощью COMSOL Multiphysics®

Инструменты моделирования COMSOL Multiphysics^® часто используются для анализа нестационарных тепловых процессов в твёрдых телах. Нестационарные модели теплообмена довольно просто настраиваются и решаются, но работа с такими моделями всё же не лишена некоторых трудностей. Скажем, интерпретация результатов может озадачить даже опытных пользователей COMSOL^®. В этой статье блога мы рассмотрим простую модель нестационарного нагрева, на примере которой мы постараемся поглубже разобраться в некоторых особенностях решения таких задач.

Простая задача о нестационарном нагреве

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

Похожая постановка задачи используется в учебной модели лазерного нагрева кремниевой пластины «Laser Heating of a Silicon Wafer», но не будем забывать, что методы и алгоритмы, обсуждаемые в этой статье, применимы и к любым другим задачам, связанным с нестационарным нагревом.

Рисунок 1. Схема расчётной области, на которой показана область нагрева.

Конечно, велик соблазн сразу построить точную 3D модель, показанную на рисунке 1, но мы начнём с более простой задачи. Как видно на рисунке 1, геометрическая модель и граничные условия в задаче обладают осевой симметрией, поэтому логично предположить, что решение также будет осесимметричным. Это означает, что модель можно упростить за счёт перехода к двумерной осесимметричной постановке задачи. (О том, как условия симметрии позволяют уменьшить размер расчётной области, можно прочесть здесь.)

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

Рисунок 2. Двумерная осесимметричная модель, эквивалентная трёхмерной модели. Показана автоматически сгенерированная расчётная сетка.

Мы рассматриваем случай, когда плотность теплового потока изменяется мгновенно; при t = 0.25 с величина теплового потока резко уменьшается. Для задания таких ступенчатых изменений потока следует использовать специальный интерфейс Events, как описано в соответствующей статье базы знаний COMSOL о решении задач с нестационарными граничными условиями. Короче говоря, интерфейс Events позволяет указать, когда именно происходит изменение теплового потока, благодаря чему решатель автоматически корректирует временной шаг в соответствующие моменты времени. Также мы хотели бы знать величину временных шагов, которые использует решатель, чтобы настроить вывод результатов на каждом шаге по времени и построить график изменения температуры в центре верхней грани цилиндра, как показано на рисунке 3.

Рисунок 3. График изменения температуры в точке показывает, что при резком изменении теплового потока решатель уменьшает шаг по времени.

Теперь мы выполним расчёт несколько раз при разных значениях относительной точности решателя (относительная точность или относительный допуск определяет критерий сходимости решателя) и сравним полученные результаты на графике (рис. 4). Мы видим, что при уменьшении относительного допуска решение ожидаемо быстро сходится к одним и тем же значениям.

Рисунок 4. График изменения температуры во времени в точке для трёх значений относительного допуска решателя.

Еще одна величина, которую можно рассчитать, — это общее количество энергии, поступающей в расчётную область. Мы проинтегрируем плотность теплового потока ht.nteflux по внешним границам и с помощью оператора timeint() рассчитаем интеграл по времени, который и даст нам изменение полной энергии. В таблице ниже приведены результаты интегрирования при трёх значениях относительного допуска нестационарного решателя. (Совет: подробности об интегрировании по пространству и времени можно узнать в этой статье базы знаний COMSOL, а о том, как рассчитать баланс энергии мы рассказывали в статье блога «Как проверить выполнение баланса массы и энергии»).

Эти данные показывают, что точность расчёта полной энергии тела практически не зависит от относительного допуска нестационарного решателя. На первый взгляд, это говорит о фантастической точности нашей модели. Однако важно отметить, что этот «эффект» является фундаментальным математическим свойством метода конечных элементов (МКЭ). Проще говоря, баланс полной энергии всегда будет соблюдаться с очень хорошей точностью. Это не означает, что расчётная модель не имеет погрешности, просто источники этой погрешности проявляются по-разному… так что давайте их поищем.

Ошибки: легко сделать, трудно оценить

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

Ошибки ввода

Ошибки ввода, как следует из названия, — это ошибки, допущенные при задании исходных параметров задачи, например свойств материала или геометрических размеров. Одной из самых серьёзных ошибок ввода является пропуск необходимых данных, скажем, когда не задано необходимое граничное условие. Не следует путать ошибки ввода с погрешностью исходных данных, которая обусловлена, например, отсутствием точных сведений о свойствах материала. Ошибки ввода можно устранить, внимательно и скрупулёзно настраивая расчётную модель, тогда как погрешность исходных данных можно оценить с помощью специального модуля «Оценка неопределённости». Будем считать, что в нашей модели ошибки ввода и погрешность исходных данных отсутствуют.

Погрешность геометрической дискретизации

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

Погрешность дискретизации решения

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

Погрешность интегрирования по времени

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

Погрешность интерпретации результатов

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

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

Итак, определившись с терминами, можем теперь вернуться к нашей модели.

Оценка погрешности в пространстве и времени

Выше мы рассматривали изменение температуры в одной точке расчётной области и обнаружили, что решение очень хорошо сходится при уменьшении величины относительного допуска нестационарного решателя. Таким образом, мы можем сделать вывод, что чем меньше относительной допуск нестационарного решателя, тем меньше погрешность решения на каждом шаге по времени. Теперь давайте проанализируем найденное пространственное распределение температуры. Для начала рассмотрим изменение температуры вдоль осевой линии. Построим график температуры для самого большого относительного допуска 0.01 в начальный момент времени и после выполнения первого шага по времени.

Рисунок 5. График изменения температуры вдоль осевой линии в начальный момент времени и на первом шаге по времени.

Как мы видим на графике, температура вдоль осевой линии не совпадает с заданным начальным распределением — в некоторых точках температура оказывается ниже начального значения. Это связано с тем, что при решении нестационарных задач в COMSOL Multiphysics используется так называемая согласованная инициализация consistent initialization. Это процедура, которая корректирует начальное распределение искомой величины так, чтобы оно соответствовало заданным начальным и граничным условиям задачи. Фактически, в рамках этой процедуры выполняется расчёт для очень маленького временного шага, который можно интерпретировать как начальный момент времени. Согласованную инициализацию можно отключить в окне настройки нестационарного решателя, а также в окне настройки узлов Explicit Events и Implicit Events, но делать это следует с осторожностью. В общем случае в мультифизических задачах, особенно в связанных с моделированием гидродинамики, использование этой процедуры может давать более надёжные результаты, поэтому в нашей модели мы не будем отключать эту опцию.

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

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

Для краткости мы пропустим вывод аналитического решения и просто приведём результат: когда теплота подводится к поверхности тела, температура на поверхности начинает расти, и внутренняя область также нагревается. Заметим, что для нагрева точек, расположенных дальше от границы, требуется больше времени. Температура внутри пластины будет изменяться неоднородно: температура в точках, расположенных дальше от поверхности, изменяется меньше, чем в точках ближе к поверхности. Важно отметить, что пространственная неоднородность температуры со временем сгладится благодаря диффузионной природе теплопроводности. Теперь давайте вернемся к нашей модели и посмотрим, как повысить её точность.

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

Рисунок 6. В последовательность команд построения сетки добавлена команда генерации погранслойных элементов на одной из граней верхней поверхности пластины.

Теперь повторно запустим расчёт, а после его завершения построим график распределения температуры в начальный момент времени и на первом временном шаге.

Рисунок 7. График распределения температуры вдоль осевой линии на начальном и первом временных шагах при использовании погранслойной сетки.

На рисунке 7 можно заметить, что колебания температуры в начальный момент времени теперь локализованы в небольшой области пространстве. Оказывается, что использование более плотной сетки приводит также к уменьшению временных шагов нестационарного решателя. Таким образом, сгустив сетку, мы уменьшили как погрешность пространственной дискретизации, так и погрешность интегрирования по времени.

Давайте посмотрим на распределение температуры по верхней границе расчётной области. На рисунке 8 показаны графики для начального момента времени и первого временного шага, полученные при относительном допуске решателя 0.01. На этих графиках мы видим довольно резкие колебания температуры в пространстве. Это сигнализирует о проблеме пространственной дискретизации. Напомним, что плотность теплового потока меняется по радиусу ступенчато, и эффект, который мы наблюдаем, отчасти напоминает эффект Гиббса.

Рисунок 8. График распределения температуры по верхней поверхности в начальный момент времени и на первом временном шаге, полученного на погранслойной сетке.

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

Рисунок 9. Окно настройки узла Mesh и расчётная сетка со сгущением в точке ступенчатого изменения плотности теплового потока.

Как видно на рисунке 10, колебания температуры уменьшились и не так сильно распространяются в пространстве и времени. Даже при относительном допуске решателя 0.01 точность решения значительно повысилась.

Рисунок 10. После сгущения сетки и задания относительного допуска 0.01 точность расчёта температуры на верхней поверхности в начальный момент времени и на первом временном шаге повысилась.

Можно можем продолжить снижать относительный допуск решателя и сгущать расчётную сетку, однако уже сейчас видно, что погрешность быстро снижается — даже оставшиеся колебания решения затухают в пространстве и времени благодаря диффузионной природе нестационарной теплопроводности. Действительно, теперь стоит затратить те же усилия, чтобы оценить погрешность исходных данных модели и её влияние на результаты моделирования.

Что ещё нужно учесть?

В этом примере заданная зона нагрева не перемещается по верхней поверхности, поэтому использованный нами приём с разбиением верхней границы на две области является разумным. Если бы расположение зоны нагрева менялось во времени, тогда пришлось бы сгустить сетка на всей верхней грани. (Прочесть о трёх методах моделирования подвижных граничных условий в COMSOL^® можно здесь.)

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

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

Заключительные замечания

В этой статье блога мы рассмотрели типовую задачу моделирования нестационарного теплообмена. Мы отметили, что существуют определённые источники погрешности, приводящие к возникновению ошибок в решении при резких изменениях теплового потока как в пространстве, так и во времени. Теперь читатели должны понимать, что это за ошибки, и знать, что они являются принципиальной особенностью метода конечных элементов, который, как и все численные методы, позволяет получить лишь приближённое решение. Хотя эти ошибки могут казаться значительными, их величина уменьшается в пространстве и времени из-за диффузионной природы уравнения теплопроводности. Мы показали, что сгущение сетки уменьшает погрешность пространственной дискретизации, а также косвенно снижает погрешность интегрирования по времени. Наконец, мы обсудили, как дополнительно сократить погрешность интегрирования по времени с помощью уменьшения относительного допуска решателя.

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

Теперь мы сможем избежать ошибок интерпретации результатов. Это позволит нам уверенно и легко строить сложные модели на основе более простых задач.

Дальнейшие шаги

Попробуйте поработать с описанной в этой статье моделью самостоятельно, кликнув по кнопке, которая приведет вас в галерею приложений: