CC BY
23
УДК 004.422
DOI: 10.25206/1813-8225-2018-160-169-173
с. В. ФЕДОРОВ
Омский государственный технический университет, г. Омск
СНИЖЕНИЕ ПОГРЕШНОСТЕЙ ПРИ РЕШЕНИИ ЗАДАЧ ТЕПЛОПРОВОДНОСТИ МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ
В статье представлен метод выбора шага по времени при решении нестационарных задач теплопроводности. Произведено сравнение с аналитическим решением. Предложен алгоритм выбора шага в автоматическом режиме. Применение данной методики позволяет существенно сократить затраты машинного времени на проведение расчета нестационарного температурного поля без ущерба в точности.
Указанная методика реализована в программе «Temper-3d», используемой для выполнения теплотехнических расчетов ограждающих конструкций.
Ключевые слова: метод конечных элементов (МКЭ), теплотехнический расчет МКЭ, «Temper-3d», валидация EN ISO 1021, автоматическая генерация конечно-элементной сети.
Введение. Расчет температурных полей методом конечных элементов широко используется в машиностроении, строительстве автомобильных дорог, промышленных и гражданских зданий, строительстве холодильных установок и др.
Этот метод также используется при производстве работ в геологоразведке, нефтедобыче, а также при расчете конструкций на пожаробезопасность.
Расчет температурных полей аналитическим методом возможно только для тел простой формы [1, 2].
Однако сложность проведения расчета, и прежде всего численная реализация, не позволяет производить расчеты реальных трехмерных конструкций за приемлемое время с требуемой точностью.
Доказано, что при измельчении конечно-элементной сети точность расчета увеличивается [3], однако это ведет к существенному увеличению размерности системы линейных уравнений, что, в свою очередь, значительно увеличивает затраты машинного времени на проведение расчета.
Размерность конечно-элементной сети для реальных трехмерных расчетов может достигать 1 млн узлов и более, а затраты машинного времени могут составлять несколько часов.
Чрезмерное уменьшение конечно-элементной сети не позволяет детально идеализировать расчетную область, а самое главное — вносит существенную погрешность.
Степень дискретизации в значительной степени зависит от вида конечного элемента, вернее, от закона распределения температуры в конечном элементе.
В данной работе используется восьмиузловой изопараметрический элемент с трилинейным законом распределения температуры. Использование данного закона распределения означает, что если температурное поле будет незначительно отличаться от линейного, то погрешность расчета будет
невелика, в случае существенно нелинейного распределения температуры необходима более мелкая дискретизация.
Например, при расчете трехслойной ограждающей конструкции достаточно всего трех конечных элементов для получения точного решения, так как распределение температурного поля в данном случае будет линейным, дополнительное измельчение в данном случае совершенно не требуется.
В случае же расчета области в виде угла требуется значительное измельчение в зоне угла, так как в этой области температурное поле существенно нелинейно.
Поэтому актуальным является вопрос об оценке погрешностей вносимых в расчет в зависимости от степени дискретизации. Также важным является разработка методики, позволяющая производить автоматическую генерацию конечно-элементной сети, что значительно сокращает затраты на подготовку исходных данных и снижает ошибки при вводе исходных данных, так как человеческий фактор сведен к минимуму.
1. Пример автоматической генерации конечно-элементной сети программой Temper-3d. Для тестирования алгоритма необходимо точное (аналитическое) решение, но таких решений найти в литературе не представилось возможным, так как решение приводилось для тел простой формы со сравнительно однородными граничными условиями. А это не представляло интереса с точки зрения метода конечных элементов, так как требовалось произвести сравнение на расчетной области с граничными условиями, которые формировали бы существенно нелинейное температурное поле.
Поэтому в качестве аналитического решения была выбрана задача, в которой в силу симметрии можно было бы произвести сравнительный анализ.
Расчетная область представляет из себя пластину размером 500x500 мм. Она представлена на рис. 1.
Рис. 1. Расчетная область
Рис. 3. Температурное поле с равномерной конечно-элементной сетью
Рис. 2. Температурное поле
Рис. 4. Расчетная область в режиме «Искажения»
Произведем расчет температурного поля методом конечных элементов по программе Temper-3d [4, 5] плоская задача.
Расчет производился при следующих граничных условиях:
— теплопроводность пластины ^ = 0,01 Вт/(мК);
— коэффициент теплоотдачи а =100 Вт/(м2К);
— температура в зоне подвода тепла (сверху и снизу) ^ = 100 °С;
— температура в зоне отвода тепла (слева и справа) ^= — 100 °С.
На рис. 2 приведено температурное поле.
Как видно из рис. 2, температурное поле симметрично, на диагоналях температура равна 0 °С.
Максимальная температура £ =99,93 °С наблю-
А * тах
дается в середине области сверху и снизу. Минимальная температура £тП =—99,93 °С наблюдается в середине области слева и справа.
Значения максимальных и минимальных температур не анализировались.
Температура во всех четырех углах равна 0 °С.
Для того, чтобы выяснить, как необходимо разбивать расчетную область на конечные элементы, для получения корректного результата, произведем расчет с равномерной конечно-элементной сетью с шагом 12,5 мм, всего получится 1600 КЭ (рис. 3).
Как видно из рис. 3, температурное поле в углах конструкции некорректно, кроме того, в точках,
близких к угловым, наблюдаются некорректные значения температур 124,66 °С (сверху и снизу), так и —124,66 °С (справа и слева), что, конечно же, может быть приемлемым, так как значения этих температур находятся за пределами температур окружающей среды 100,00 °С, -100,00 °С.
Причем температурное поле в центре области отображено вполне корректно.
Можно предположить, что более мелкая дискретизация в угловых зонах приведет к более корректному расчету, т.е. использовать неравномерную конечно-элементную сеть.
Однако априори сложно определить, как необходимо производить дискретизацию, чтобы получить корректное решение.
Для решения этой проблемы в программе Temper-3d, используется понятие «искажение» в КЭ, которые рассчитываются для каждой грани конечного элемента, эти искажения представляют из себя отклонение температур для четырех узлов каждой грани конечного элемента от плоскости [6].
На рис. 4 приведена расчетная область, где изображены искажения.
Как видно из рис. 4 наибольшие искажения возникают вблизи углов и достигают значения 7,67 °С.
Как показали результаты численных исследований, для получения корректного расчета необходимо, чтобы искажения не превышали 1 °С.
Рис. 5. Температурное поле при использовании неравномерной конечно-элементной сети
Программа производит автоматическую дискретизацию области на конечные элементы таким образом, чтобы в любой грани конечного элемента искажения не превышали наперед заданного значения, например 1 °С.
На рис. 5 приведено температурное поле, полученное при автоматическом генерировании конечно-элементной сети (представлен фрагмент расчетной области, всего 4620 КЭ).
Как видно из рис. 5, размер конечного элемента в центре конструкции достаточно велик (25 мм), а по мере приближения к углам области размер зна-
чительно уменьшается и составляет в углу 0,2 мм, отношение максимального конечного элемента к минимальному в данном расчете составляет 125.
Максимальные искажения для данного расчета составляют 2,25 °С, меньшее значение получить не удалось, так как в программе стоит ограничение на минимальный размер КЭ (0,15 мм).
Данное обстоятельство связано со специфичностью расчета (в углу сверху температура 100 °С, а слева -100 °С).
При проведении реальных расчетов удается получить практически сколь угодно маленькое значение искажения, вплоть до 0,01 °С.
2. Возможности программы Temper-3d в соответствии с требованиями EN ISO 10211. Программа Temper-3d предназначена для выполнения теплотехнических расчетов ограждающих конструкций с применением температурных полей и тепловых потоков через узлы ограждающей конструкции. При создании плоских и объемных моделей узлов ограждающей конструкции в Temper-3d производится автоматическое разбиение конечно-разностной сетки, соответствующей требованиям EN ISO 10211 [4, 5].
Для проверки точности расчета имеется возможность производить измельчение сетки в кратное количество раз для оценки корректности результатов согласно требованием EN ISO 10211. Ниже приведен пример расчета плоской ограждающей конструкции.
Всего проведено девять расчетов, из них — восемь с равномерной конечно-элементной сетью,
Рис. 6. Дискретизация на конечные элементы (расчеты № 1-№ 3)
Рис. 7. а — расчетная область, расчет № 4; б — температурное поле, расчет № 9
(7 Измельчать автоматически Параметр измельчения
rarer
1000
- Искажения -
г г с г г Fs
Измельчить для проверки
Рис. 8. Панель настройки для проведения расчетов
один — с неравномерной конечно-элементной сетью при автоматическом разбиении. Во всех расчетах использованы идентичные граничные условия: температура наружного воздуха —30 °С, температура внутреннего воздуха 20 °С.
Расчетная область, для расчета № 1 состояла из пяти конечных элементов, каждый размером 640x640 мм, для расчета № 2 из 20 конечных элементов каждый размером 320x320 мм, и т.д. Каждый последующий расчет получался путем увеличения в 2 раза в направлении х и у количества конечных элементов.
На рис. 6 показана дискретизация (разбиение) на конечные элементы для трех расчетов.
На рис. 7а приведена дискретизация (разбиение) расчетной области на конечные элементы для расчета № 4 с указанием размеров модели. Рисунки для расчетов № 5 — № 8 не приведены.
В качестве конечного элемента использовался восьмиузловой изопараметрический конечный элемент [4, 5]. Расчеты производились на удаленном сервере по сети Интернет [7].
На рис. 7б приведено температурное поле для расчета № 9 (промежуточный расчет с неравномерной конечно-элементной сетью). В качестве контрольного значения принята температура внутреннего угла в точке А.
Результаты расчета представлены в табл. 1.
Как видно из табл. 1, абсолютная погрешность расчета контрольной температуры уменьшается с увеличением количества конечных элементов, что свидетельствует о корректности проведенных расчетов. Также уменьшается значение Af в расчете № 4 (AfRs. = 0,005A). Значит, этот и последующие расчеты проведены корректно (Условие 1 по EN ISO 10211) [8, 9].
Погрешность тепловых потоков при измельчении в 2 раза уже в расчете 3 менее 1 %. Это значит, что все последующие расчеты после № 3 корректные (Условие 2 по EN ISO 10211).
Погрешность входящих и выходящих тепловых потоков в табл. 1 не приведена, так как для всех расчетов она значительно меньше 0,01 % (Условие 3 по EN ISO 10211).
Корректным также является расчет с неравномерной конечно-элементной сетью при меньшем количестве конечных элементов.
При решении данной задачи корректным следует считать расчет № 4 и все последующие, так как именно в этих расчетах выполняются все три условия.
Таблица 1
№ расчета Количество КЭ в области 640x640 мм, лхл Максимальное искажение температуры, ° С Температура в точке А Температурный коэффициент fRsi в точке А Изменение температурного коэффициента Л fRi Количество КЭ Входящий тепловой поток Q, Вт Приведенное сопротивление теплопередачи R0 (м2х°С)/Вт Изменение теплового dQ, % Погрешность температуры в точке А, dt, °С*
1 1 27,96 8,89 0,778 5 14,41 0,888 0,81
2 2 11,07 9,01 0,78 0,002 20 14,26 0,897 1,0057 0,93
3 4 4,96 8,49 0,77 -0,01 80 14,2 0,901 0,4196 0,41
4 8 2,3 8,24 0,765 -0,005 320 14,18 0,903 0,1429 0,16
5 16 1,06 8,14 0,763 -0,002 1280 14,18 0,903 0,0434 0,06
6 32 0,47 8,1 0,762 -0,001 5120 14,17 0,903 0,0123 0,02
7 64 0,21 8,09 0,762 0 20480 14,17 0,903 0,0033 0,01
8 128 0,09 8,08 0,762 0 81920 14,17 0,903 0,0009 0
9 - 0,07 8,13 0,763 0 2972 14,17 0,903 -0,0027 0,05
* при определении абсолютной погрешности в качестве точного решения принят расчет № 8
Искажения температуры контрольной точки в расчете № 4 составляют 2,3 °С. Однако, как показала практика расчета трехмерной модели, результаты будут корректными при искажениях от 0,6 до 1 °С. В этой связи критерием корректности расчета не может служить величина искажений температуры. Она необходима только для генерации рационально конечно-элементной сети.
2. Рекомендации по выбору шага. В используемых в настоящее время нормативных документах, для расчета ограждающих конструкция МКЭ рассматривается вопрос о корректности расчета всей ограждающей конструкции, включая внутренние точки только с точки зрения интегральных характеристик, таких как, например, входящий и выходящий тепловые потоки [9, 10].
Вопрос о погрешности расчета температурных полей для точек, не принадлежащих внутренней поверхности ограждающей конструкции, EN ISO 10211 не регламентирован.
В этой связи рекомендуется следующий алгоритм проведения расчета температурного поля и теплового потока через узел ограждающей конструкции с помощью программы «Temper-3d»:
1. Произвести расчет с автоматической дискретизацией на конечные элементы ограждающей конструкции, при этом в программе устанавливаются искажения 1,5 °C и убедиться, что погрешность входящих и выходящих потоков не превышает 0,01 %. Данный параметр указан в файле *.nmn (Условие 3).
2. Произвести повторный расчет ограждающей конструкции, измельчив конечно-элементную сеть в 2 раза в по каждой координате. В случае трехмерного расчета количество конечных элементов увеличится в 2*2*2 = 8 раз.
3. Сравнить значения тепловых потоков и убедиться, что при измельчении сети в два раза по каждой координатной оси изменение теплового потока не превышает 1 % (Условие 1).
4. Произвести расчет температурного коэффициента внутренней поверхности, используя в качестве расчетных узлов точки на внутренней поверхности. В этом случае отличие коэффициента не должно превышать 0,005 для любой точки, принадлежащей внутренней поверхности. Например, если температура внутреннего воздуха 20 °С, а наружного -30 °С, то разница температур в точках на внутренней поверхности не должна превышать 0,25 °C при измельчении конечно-элементной сети в два раза по каждой координатной оси (Условие 2).
5. Если выполнены условия 1, 3 и 4, то результатами расчета можно воспользоваться, в противном случае необходима более мелкая дискретизация, для этого надо повторить все с пункта №1, при этом установив меньшее значение параметра искажения, например, 1 °C и так далее, пока не будут выполнены все три условия корректности. На рис. 8 приведена панель настройки для расчетов. С ее помощью можно установить требуемые искажения,
параметр измельчения и произвести дополнительное измельчение конечно-элементной сети для проведения контрольного расчета.
6. Считать расчет корректным, если искажения менее, например, 1 °C; не соответствует нормативным требованиям EN ISO 10211, этот параметр используется только для автоматической генерации рациональной конечно-элементной сети.
Библиографический список
1. Лыков А. В. Теория теплопроводности. М.: Высшая школа, 1967. 600 с.
2. Лыков А. В. Тепломассообмен. 2-е изд., перераб. и доп. М.: Энергия, 1978. 480 с.
3. Зенкевич О. Метод конечных элементов в технике: мо-ногр. М.: Мир, 1975. 538 с.
4. Теплотехнические расчеты в «Temper-3d». URL: http:// www.temper3d.ru (дата обращения: 01.02.2018).
5. Федоров С. В. Свидетельство об официальной регистрации программ «Temper-3d» для ЭВМ № 2006610359 // Реестр программ для ЭВМ. М.: ФИПС. № 2005613080 от 20.01.2006.
6. Федоров С. В. Алгоритм автоматической генерации конечно-элементной сети для теплотехнических расчетов // Информационный бюллетень Омского научно-образовательного центра ОмГТУ и ИМ СО РАН в области математики и информатики. 2017. Т. 1, № 1. С. 128-129.
7. Федоров С. В. Клиент-серверные расчеты по e-mail // Прикладная математика и фундаментальная информатика. 2016. № 3. С. 204-207.
8. Федоров С. В., Терехова И. А. оценка корректности теплотехнических расчетов ограждающих конструкций методом конечных элементов // Прикладная математика и фундаментальная информатика. 2017. Т. 4, № 1. С. 31-42.
9. EN ISO 10211:2007 Тепловые мостики в зданиях. Тепловые потоки и температура поверхности. Подробные расчеты. Авеню Marnix 17, B-1000 Брюссель, 2009. 64 с.
10. СП 50.13330.2012 Тепловая защита зданий (взамен СНиП 23-02-2003). Введ. 2013-07-01. М.: Госстрой России, 2013. URL: http://docs.cntd.ru/document/1200095525 (дата обращения: 01.02.2018).
ФЕДОРОВ Сергей Витальевич, кандидат технических наук, доцент (Россия), доцент кафедры «Прикладная математика и фундаментальная информатика».
SPIN-код: 9085-8584, АиШогГО (РИНЦ): 903562 Адрес для переписки: temper99@mail.ru
Для цитирования
Федоров С. В. Снижения погрешностей при решении задач теплопроводности методом конечных элементов // Омский научный вестник. 2018. № 4 (160). С. 169-173. Б01: 10.25206/1813-8225-2018-160-169-173.
Статья поступила в редакцию 11.04.2018 г. © С. В. Федоров
