CC BY
55
ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ
УДК 519.6:697.1
Н. Н. Павлов, В. Ю. Шадрин
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ И АНАЛИЗ РАСЧЕТА ТЕПЛООБМЕНА В НАРУЖНЫХ ОГРАЖДАЮЩИХ КОНСТРУКЦИЯХ С ОСОБЕННОСТЯМИ
Посвящена разработке и анализу численных методов решения сеточных уравнений, аппроксимирующих третью краевую задачу для двумерного уравнения Лапласа. Показана эффективность метода сопряженных градиентов для расчета стационарных температурных полей наружных ограждающих конструкций зданий и сооружений с угловыми элементами, состоящих из материалов с тысячекратными разностями значений коэффициентов теплопроводности.
Ключевые слова: ограждающие конструкции, математическое моделирование, разностная схема, метод сопряженных градиентов, погрешность аппроксимации.
N. N. Pavlov, V Yu. Shadrin
Numerical modeling and analysis of heat exchange calculation in outer bounding constructions with peculiarities
The article is about working out and analysis of digital method of finite-difference equation solution approximating the third boundary value for two-dimensional Laplace equation. The effectiveness of the conjugate gradient method for computation of stationary temperature fields of outer bounding buildings’ constructions and constructions with angular elements consisting of materials with thousand-fold differences of heat passage values is showed.
Key words: bounding constructions, math modeling, differential scheme, the conjugate gradient method, approximation error.
Данная работа посвящена разработке и исследованию численных методов расчета двумерного стационарного температурного поля наружных ограждающих конструкций с особенностями.
Как показал опыт работы с программным продуктом «SHADDAN» [1], разработанным авторами, возникают некоторые проблемы при расчете температурных полей наружных ограждающих конструкций с угловыми элементами и наличия в ограждениях соседних материалов с несоизмеримыми размерами и
ПАВЛОВ Никифор Никитич - к. ф.-м. н., доцент кафедры информационных технологий Института математики и информатики СВФУ им. М.К. Аммосова.
E-mail: pnn10@mail.ru
ШАДРИН Василий Юрьевич - к. ф.-м. н., профессор кафедры высшей математики Института математики и информатики СВФУ им. М.К. Аммосова.
E-mail: vshadr@mail.ru
тысячекратным перепадом значений коэффициентов теплопроводности. При расчете таких ограждений необходимо увеличить количество узлов расчета, что ведет к удлинению времени расчета, увеличению требований к оперативной памяти и контролированию конфигурации сетки.
Для решения систем уравнений, возникающих при дискретизации третьей краевой задачи для двумерного уравнения Лапласа с переменными коэффициентами, в программе «8НАЭБА№> применялся итерационный метод Зейделя, который неплохо зарекомендовал себя для достаточно широкого класса наружных ограждающих конструкций и был выбран из-за простоты реализации и удобства для организации интерфейса.
В данной работе для решения систем уравнений предложен итерационный метод сопряженных градиентов [2], проведен анализ применения и сравнения с методом Зейделя для расчета температурного поля наружных ограждающих конструкций с особеннос-
тями.
1. Постановка задачи
Рассмотрим третью краевую задачу для двумерного уравнения Лапласа с кусочно-постоянными коэффициентами теплопроводности:
дТ
дТ
= 0 , (х1, х2) еО,
лдТ+аср (Г - Тср )=0, (X! , х2) єда,
(1)
(2)
где Т (х1, х2) - температура точки конструкции с координатами ( Х1,Х2 ) , Л = Л(х1,Х2) - коэффициент теплопроводности, ЭО - граница области О, а - ко-
эффициент теплообмена со средой, Тср - постоянная
дт
дп
температура внешней среды,
производная по
внешней нормали. Область О состоит из конечного числа прямоугольников со сторонами, параллельными координатным осям, как показано на рис. 1. Каждый прямоугольник имеет свой коэффициент теплопроводности Л .
Для решения задачи (1)-(2) применим метод сеток.
Введем неравномерную сетку
Ък = {(Xа , х2а), іа= 0,1,..., Иа, а = 1,2}, щ =ак + ук,
Г\
.... н П И н г- / \уЬ.
1 .
Пг
ТІ '
— — — и «—
Рис. 1. Сеточная область
где ($ - множество внутренних узлов, у - мно-
Н Н
жество граничных узлов, = х1а — х1а 1,
К = %= 0.5( Н + к1/1).
На рис. 1 знаком о обозначены внутренние узлы, а знаками х, р, Ц ^, ^, ^, ^, ^, ^ обозначены граничные узлы. Для множества граничных узлов, обозначенных знаком X, введем множества левых
и правых граничных узлов у^ 1 и у1аг), а = 1,2.
Обозначим ^<+) = уг + Уп + Уп + У г,
7<-) = У# + У" + у# + у".
В области а рассмотрим сеточные функции
у== у(\, у^ =мх(±1з-аЬ
и(±) = и( Х^3-а )).
Дискретизацию задачи проведем с помощью интегро-интерполяционного метода, который приводит к консервативной разностной схеме со вторым порядком аппроксимации по обеим пространственным переменным [3]:
Аи = / , (3)
где оператор А = А + А2, / = / + /2:
12((^<а+ % ) Х + (аа У г ) Х )> х ’
2 ха ха ха ха
2(Г(а«+)Уха -«<+)У) + ^-(а?Уха -у)), х
Ь1 (-«<+> у - «<+> ух ) + К-(-а<-> у - а? Ух )), 2 Па К
1Кг(а<‘+)Уха -а<+у')’ хе г#>
а
2т(а°}Уха)-а<-у)’х’
2К-(-«<+)У-а(а+ Ух )’ х
2 П а Х а
\Ь(-а( )У-№Ух )’ хе ^’
2 К а ха
2(К-(а*}Уха -У) + К+ Ух )х.), хе 7#,
2 К а а Ла ха
2(-Г-(аа+ Ух - «<+) У) + (аа} Ух )„)’ х еУ"’
2 Ка а Ха ха
2((а«+)Ух ) х + К~(-а< )У - а<а-Ух ))’ хе У!,
2 Ха Ха Па Ха
2((аа} УХ ) Х + Т-(-а<Л> У - ^ УХ ))’ Х еУг ’
2 Ха ха К Ха
АаУ =
/а =
О, х є сок,
2 Т {-а^ + а^ Т. х єУа, 2 к-( «£> + С%> )Тср , х єуа,
——а(+) Т х є і/+)
2 к а ср Тср ’ х ' ’
ІГ ^ . Х ЄГ-
Систему (3) будем решать методом сопряженных градиентов [2]:
Вуы =ак+1(В-тк+1А)ук + (1 -а^Ву1-1 + ак+1тк+1/, к = 1,2,...,
(4)
Ву1 = (В -Ті А) у0 +т/,
где ак+1 и Тк+1 - итерационные параметры. Вычисления по (4) основаны на представлении
у*+1 = ак+і ук +(1 -ак+1) у к-1 + ак+1тк+1ак.
Здесь ($к = В ~ггк.
Итерационные параметры рассчитываются по формулам:
(ок, гк) , П1
Тк+і = Л к , к = 0Д,...,
(АО,О)!
ак+1 = (1
(°, гк) 1_)-^ к = 0,1,2,...,
тк (О 1, гк 1) ак
а1 =1.
В качестве переобуславливателя В используется диагональная часть оператора А. Такой выбор особенно важен для данной задачи, т. к. коэффициент теплопроводности является сильно переменным.
2. Численные эксперименты
Проведены численные эксперименты для четырех видов наиболее распространенных типов наружных ограждающих конструкций:
2.1. Однородная деревянная стена (точное решение известно).
2.2. Неоднородная стена с большими перепадами коэффициентов теплопроводности (точное решение известно).
2.3. Неоднородная стена с теплопроводным включением.
2.4. Угловое соединение.
Введем обозначения:
е =| Тк-Тк-11 - разность между соседними итерациями Тк и Тк-1 в равномерной норме, к - номер итераций, ДР=||Р„ар-Рвнутр ||- разность между
тепловыми потоками через наружную поверхность (Р ) и через внутреннюю поверхность (Рвнутр)
ограждения, 1 - время расчета, Тнар - вычисленная
температура на наружной границе ограждения, Твнутр
- вычисленная температура на внутренней границе ограждения.
При численных экспериментах для дополнительного контроля правильности расчетов проводилось вычисление разности входящих и выходящих тепловых потоков Др через ограждение. Такой контроль особенно важен при отсутствии точного решения температурного поля. В методе Зейделя увеличение количества узлов сетки ведет к возрастанию разности входящих и выходящих тепловых потоков через граничные поверхности ограждающих конструкций. Это является естественным следствием плохой обусловленности матрицы системы и накопления вычислительных погрешностей с возрастанием числа узлов сетки при численном интегрировании по контуру ограждения. Тем не менее, как будет видно из таблиц результатов расчетов, тепловой баланс ограждающих конструкций выполняется с достаточной точностью.
2.1. Рассматривалась прямая вертикальная деревянная стена толщиной 0,2 метра, длиной 1 метр и коэффициентом теплопроводности Х=0,14 Вт/(м^°С). Точное решение легко вычисляется, температуры на наружной и внутренней границах соответственно
равны Т =-51,945 °С, Т =15,568 °С. Результаты
А нар 7 7 внутр 7 ^
расчетов приведем в табл. 1-5.
Таблица 1
Результаты расчета методом Зейделя для е=10-3
сетка к дд і Т нар Т внутр
10Ч50 163 7.20^10-2 <1 с -51.943 15.572
20Ч100 533 3.5810-1 <1 с -51.937 15.585
40Ч200 1579 1.5110+0 3 с -51.911 15.648
100Ч500 5313 9.66^10+0 1 мин. -51.730 16.104
200Ч1000 8795 3.44^10+1 6 мин.30 с -51.078 17.222
Т
Таблица 2
Результаты расчета методом Зейделя для е=10-4
сетка к дд 1 Т нар Т внутр
10Ч50 227 7.58-10'3 <1 с -51.945 15.568
20Ч100 765 3.50-10-2 <1 с -51.944 15.570
40Ч200 2502 1.49 10-1 4 с -51.942 15.576
100Ч500 11071 9.60-10-1 2 мин. -51.924 15.621
200Ч1000 30471 3.87-10+° 23 мин. -51.860 15.786
Таблица 3
Результаты расчета методом Зейделя для е=10-5
сетка к дд 1 Т нар Т внутр
10Ч50 285 7.41-10"4 <1 с -51.945 15.568
20Ч100 996 3.45-10-3 <1 с -51.945 15.568
40Ч200 3425 1.4810-2 6 с -51.945 15.569
100Ч500 16825 9.56-10'2 3 мин. -51.943 15.573
200Ч1000 53471 3.87-104 40 мин. -51.937 15.590
Таблица 4
Результаты расчета методом Зейделя для е=10-6
Как видно из таблиц, для метода Зейделя е = 10-3 недостаточно для получения точных результатов на любой сетке. Тем не менее, он достаточно хорошо работает для е=10-6, хотя не следует увлекаться слишком мелким разбиением области расчета. Результаты табл. 5 свидетельствуют о большой скорости сходимости метода сопряженных градиентов, когда уже при е =10-3 достигается точное решение и число итераций линейно зависит от числа узлов сетки.
2.2. Рассмотрим температурное поле ограждения, которое состоит из плоского пенополистирола (Х2=0,041 Вт/(м^°С)) толщиной 10 сантиметров,
облицованного с обеих сторон алюминиевой пластиной (Х1=221 Вт/(м^°С)) толщиной 2 миллиметра.
Перепад значений коэффициентов теплопроводности составляет 1/^=221/0,041=5390,24. Точное решение также легко вычисляется, температуры на наружной и внутренней границах соответственно равны Т =-52,745 °С, Т =17,681°С. Сетка фиксирована:
нар 7 7 внутр 7 'гг
104^10, то есть шаг сетки по оси Ох И =0,001, шаг
7 х 7 7
сетки по оси Оу Ь=0,01.
Результаты расчетов приведены в табл. 6-7.
Из табл. 6 и табл. 7 видно, что метод сопряженных градиентов эффективнее метода Зейделя для случая больших градиентов температуры даже для небольшого количества узлов сетки. Очевидно, что это преимущество многократно возрастет при увеличении количества узлов.
Таблица 6
Результаты расчета методом Зейделя
сетка к дд 1 Т нар Т внутр
10Ч50 343 7.25-10-5 <1 с -51.945 15.568
20Ч100 1227 3.4110-4 1 с -51.945 15.568
40Ч200 4347 1.46 10-3 8 с -51.945 15.568
100Ч500 22577 9.5310-3 4 мин. 10 с -51.945 15.568
200Ч1000 76464 3.86-10'2 56 мин. -51.944 15.570
е к дд 1 Т нар Т внутр
10-3 19451 3.36-10+1 5 с -33.822 6.288
10-5 171502 4.46-10-1 46 с -52.752 17.189
10-7 398143 4.62-10-3 1 мин. 45 с -52.745 17.676
10-8 511490 4.62-10-4 2 мин. 15 с -52.745 17.681
Таблица 5
Результаты расчета методом сопряженных градиентов для е=10-3, 10-4, 10-5, 10-6
Таблица 7
Результаты расчета методом сопряженных градиентов
сетка к дд 1 Т нар Т внутр
10Ч50 12 8.86-10-12 <1 с -51.945 15.568
20Ч100 23 1,26'10-12 <1 с -51.945 15.568
40Ч200 43 1.06-10-10 <1 с -51.945 15.568
100Ч500 103 1.59-10-10 2 с -51.945 15.568
200Ч1000 203 1.83 10-9 15 с -51.945 15.568
е к дд 1 Т нар Т внутр
10-3 106 5.5110-6 <1 с -52.745 17.681
10-5 106 5.5110-6 <1 с -52.745 17.681
10-7 106 5.5110-6 <1 с -52.745 17.681
10-8 106 5.5110-6 <1 с -52.745 17.681
*1 ^ ^1 0,002 1 ! 0,20
^3 0,19 1
/ ! 0,01
К ТИ(=21°С
Рис. 2. Ограждение с теплопроводным включением
2.3. Теперь рассмотрим случай расчета двумерного температурного поля ограждения с теплопроводным включением, изображенного на рис. 2.
В горизонтальном слое из пенополистирола (Х2=0,041 Вт/(м^°С)) имеется тонкий вертикальный металлический стержень (Х1=58 Вт/(м^°С)). Точное решение неизвестно, о результатах расчета можно судить лишь по качеству температурного поля. Размеры ограждения указаны на рис. 2. Сетка в расчетах бралась постоянной: 42*40.
Введем обозначения:
Т - максимальная температура наружной по-
макс.нар. 1 ^ 1 1 ^
верхности, которая, очевидно, находится посередине наружной поверхности, Т - минимальная темпе-
1 ^ 1 у мин.вн.
ратура внутренней поверхности, которая, очевидно, находится посередине внутренней поверхности ограждения.
В этом более сложном случае, хотя число итераций несколько возросло для меньшего, чем в предыдущем примере, числа узлов сетки, метод сопряженных градиентов намного предпочтителен.
2.4. Наконец, рассмотрим наиболее сложный случай - расчет температурного поля угловой наружной ограждающей конструкции. Как известно, при проектировании зданий и сооружений важное значение придается точному прогнозированию температуры во внутренних углах помещений, так как в них чаще всего может выпадать конденсат.
Приведем результаты расчета для прямоугольного двумерного угла из сосны (X = 0,14 Вт/(м^°С)) толщиной 0,2 м и длиной 1 м по наружной поверхности. Данный пример интересен тем, что расчетная область имеет внутренний угол больший, чем 180 градусов.
Введем обозначения:
Т - вычисленная температура внутреннего угла ограждения. Сетка равномерная и состоит из NxN узлов.
Результаты расчетов приведены в табл. 10-12.
По результатам экспериментов видно, что расчет угловых соединений представляет определенные трудности. Это объясняется большими вторыми производными вблизи внутренних угловых узлов [4], которые в данном случае являются величинами порядка 3ДТ0+5. Поэтому, если требуется большая точность температур внутренних угловых точек, то необходимо внимательно следить за размерами сетки, точностью соседних итераций, разностью тепловых потоков через внутреннюю и наружную поверхности ограждающих конструкций. Естественно, что время расчета, погрешность приближенного решения будут возрастать с усложнением геометрии
Таблица 8
Результаты расчета методом Зейделя
є к Ґ Т мин.вн. Т макс.нар.
10-3 3608 8.37'10+0 2 с 18.447 -7.305
10-5 122598 2.28-10-1 49 с 17.642 -20.845
10-7 272234 2.27-10-3 1 мин. 50 с 17.615 -21.144
10-8 347052 2.27-10-4 2 мин. 19 с 17.615 -21.147
Таблица 9
Результаты расчета методом сопряженных градиентов
є к Ґ Т мин.вн. Т макс.нар.
10-3 502 2.42-10-2 <1 с 17.613 -21.147
10-5 560 4.1410-6 1 с 17.615 -21.147
10-7 608 3.72-10-8 1 с 17.615 -21.147
10-8 631 1.77-10-8 1 с 17.615 -21.147
Таблица 10
Сравнение результатов расчета для е=10-3
Метод Зейделя Метод сопряженных градиентов
N к дд X Т угл к дд X Т угл
50 192 1.610-1 <1 с 8.306 63 2.3^10-3 <1 с 8.298
100 624 7.110-1 1 с 8.002 107 3.310-3 1 с 7.967
500 7204 9.810+0 6 мин. 8.721 428 1.7^10-2 24 с 7.770
1000 14654 3.9-10+1 48 мин. 11.330 791 3.5^10-2 3 мин. 7.765
2000 19013 1.110+2 4 часа 20 мин. 13.912 1397 6.7^10-2 20 мин. 7.779
4000 - - - - 2279 9.310-2 2 часа 8 мин. 7.772
Таблица 11
Сравнение результатов расчета для е=10-6
Метод Зейделя Метод сопряженных градиентов
N к дд X Т угл к дд X Т угл
50 375 1.810-4 1 с 8.298 101 4.2^10-7 <1 с 8.298
100 1358 6.4^10-4 3 с 7.966 173 8.110"7 1 с 7.966
500 25516 8.610-3 21 мин. 7.767 702 1.6 10-5 39 с 7.766
1000 87404 3.5^10-2 4 часа 45 мин. 7.759 1369 1.7^10-5 5 мин. 7.755
2000 290847 1.510-1 75 часов 7.767 2661 2.3^10-5 35 мин. 7.751
4000 - - - - 5114 4.110-5 4 часа 32 мин. 7.750
Таблица 12
Сравнение результатов расчета для е=10-8
Метод Зейделя Метод сопряженных градиентов
N к дд X Т угл К дд X Т угл
50 496 6.810-7 1 с 8.298 116 7.4^10-9 <1 с 8.298
100 1845 2.9^10-6 4 с 7.966 197 1.610-8 1 с 7.966
500 37715 8.0^10-5 31 мин. 7.766 853 1.5^10-8 48 с 7.766
1000 136227 3.310-4 7 часов 30 мин. 7.755 1665 6.510-8 6 мин. 7.755
2000 - - - - 3276 3.110-7 47 мин. 7.751
4000 - - - - 6412 5.410-7 5 час 43 мин. 7.750
области расчета, наличием неоднородных материалов.
Данная работа нашла применение в совершенствовании программы «8НАВБА^2Б». Новая версия программы существенно улучшена не только в отношении совершенствования математического аппарата и методов решения, но и в отношении удобства пользования программы с учетом пожеланий и замечаний строителей-проектировщиков.
Л и т е р а т у р а
1. Данилов Н. Д., Шадрин В. Ю., Павлов Н. Н.
Прогнозирование температурного режима угловых соединений наружных ограждающих конструкций // Промышленное и гражданское строительство. - 2010. - № 4. С. 20-22.
2. Самарский А. А., Вабищевич П. Н. Численные методы решения задач конвекции - диффузии. - М.: Эдито-риал УРСС, 1999. - 248 с.
3. Самарский А. А. Теория разностных схем. - М.: Наука. Гл. ред. физ.-мат. лит., 1977. - 656 с.
4. Шайдуров В. В. Многосеточные методы конечных элементов. - М.: Наука. Гл. ред. физ.-мат. лит., 1989. -288 с.
