Моделирование нестационарного процесса сопряженного теплообмена горного массива и рудничного воздуха с применением технологии вычисления на графических процессорах Nvidia CUDA Текст научной статьи по специальности «Математика»

ВЕСТНИК ПЕРМСКОГО УНИВЕРСИТЕТА

2012 Серия: Физика Вып. 4 (22)

УДК 004.94+004.942

Моделирование нестационарного процесса сопряженного теплообмена горного массива и рудничного воздуха с применением технологии вычисления на графических процессорах NVIDIA CUDA

А. Р. Куцев

Пермский государственный национальный исследовательский университет, 614990, Пермь, ул. Букирева, 15

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

Ключевые слова: нестационарный процесс, №УШ1А СЦЮА, сопряженный теплообмен, горная выработка, рудничный воздух.

1. Введение

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

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

© Куцев А. Р., 2012

дения) массива, которая будет тем больше, чем дольше длится процесс прогревания (рис. 1).

Тг

Я

Рис. 1. Цилиндрическая расчетная область для отдельной горной выработки

Помимо теплопереноса в рудничном воздухе мы должны рассматривать теплоперенос в горном массиве. В данном случае существенно не только

распределение температуры вдоль выработки, но и радиальное распределение температуры в глубь горного массива.

Задача о теплопереносе интересует нас прежде всего в сетевой постановке, когда температурное поле моделируется не в отдельной выработке, а во множестве выработок с учетом согласования температур и тепловых потоков в разветвлениях выработок.

Пусть по цилиндрическому каналу радиуса Я] движется однородный однонаправленный поток газа (воздуха) с постоянной скоростью и. Толщина стенки канала Я2 - Я], причем Я2 >> Я] (рис. 2). Температура газа в начальный момент времени т=0 во всем канале равна Т], температура стенки в начальный момент времени составляет Т2, причем Т2ФТ]. Температура на внешней поверхности стенки канала в любой момент времени т>0 остается постоянной и равной Т2.

Рис. 2. Схема расчетной области

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

Так как расчетная область обладает осевой симметрией, задачу будем решать в цилиндрических координатах. Ось Oz цилиндрической системы координат направим по оси канала (рис. 2). Поперечными течениями в газе пренебрегаем, тогда, поскольку Т = Т(г,т), уравнения, описывающие изменение температуры в газе и в стенке, примут вид [1]

дТ

дт

г

■ = а

\

д 2Т 1 дТ

—2 +------------

дг г дг

\

у

(1.1)

В уравнении (1.1) коэффициент температуропроводности при г < Я1

а = -

Срг рг

а при Я1 < г < Я2

а = -

Сс Рс

(1.2)

(1.3)

Здесь Хг, с и рг — коэффициент теплопроводности, изобарная удельная теплоемкость и плотность газа, а Яс, сс и - коэффициент теплопроводности, удельная теплоемкость и плотность материала стенки.

Уравнение (1.1) дополним начальными условиями

= Т2,

Я, <г<Я-, 2 -

(1.4)

граничным условием на оси канала (условие симметрии)

дТ_

дг

= 0,

(1.5)

граничным условием на внутренней поверхности стенки канала

аг (Т - Т)

(1.6)

где Тг и Тс — температуры газа и стенки, и граничным условием на внешней поверхности стенки канала

Т = Т2

I г=R, 2

(1.7)

Граничные условия (1.5) — (1.7) должны выполняться в любой момент времени т .

Для численного решения уравнения (1.1) построим сетку с шагом А г1 по пространственной координате г в области 0 < г < Я1, с шагом Аг2 по пространственной координате г в области Я1 < г < Я2 и с шагом Ат по времени (рис. 3). Причем при г = Я1 создается двойная точка. Координаты узлов сетки

где Аг, =

гк = (к - 1)Аг1 , к = 1,Ы1 , Я,

(1.8)

относятся к газовой среде, а узлы

с координатами

Г = я1 + (к-Nl - 1)Аг2, к = N1 + 1,Ы1 + N2 , (1.9)

где А г =

Ъ

N2 - 1

, находятся в стенке канала.

Я

я

с

г = 0

г-Я

Т П + 1 Т п

k___к_

А т

= а

грп "У'Т’п і тп 7 тп Т'п

_к+1 ~ 2_к _к-1 + 1 _к+ 1 ~ _к-1 | (1 10)

( А г)г

2Аг

7'П + 1 грп . аАТ / п - п п \

к _ Тк + /л_. ,2 (Тк+1 - 2Тк + Тк-1 )

( Аг}2 а Ат

(1.11)

2гк А г

(тп - тп )

\т к+1 т к-1 ) ■

дТ

дт

дг

гг’П +1 ______ 'тЧХ

Т1 ~ Т1 +

2 (г; - 2т1 + г;).

Из условия симметрии Т0 = Т2И и для определения температуры в точке К = 1 получаем

грп+і , 2аАт (ггП грп\

' = ' 7X^71)■

(1.12)

Для определения температуры в точках k = N и k = N +1 запишем граничное условие (1.6) в конечно-разностном виде:

а.

(■тN

Аг

(1.13)

Рис. 3. Разностная сетка

Значение температуры Т(г,т) в точке гк в момент времени пАт,п = 0,1,2,... будем обозначать через

Гп

к , тогда конечно-разностный аналог уравнения (1) будет иметь вид

Таким образом, уравнение, позволяющее по известным значениям температуры на п -м временном слое определить температуру в узлах сетки на п +1 временном слое, будет иметь вид

Уравнение (1.11) применимо только для внутренних узлов сетки, т. е. для к = 2,Ы1 -1 и для

к = Ы, + 2, N2 + Ы2 -1. Причем при к = 2,Ы{ -1 в этом уравнении Аг = Ат1, гк и коэффициент температуропроводности а вычисляются по формулам (1.8) и (1.2), а при к = Ы, + 2,Ы, + Ы2 -1 в уравнении (1.4) Аг = Аг2, гк вычисляется по формуле (1.9), а коэффициент температуропроводности - по формуле (1.3).

При К = 1 на оси симметрии, с учетом граничного условия (5), уравнение (1) принимает вид

Если за пределами расчетной области ввести на расстоянии Аг1 фиктивную «нулевую» точку, то получим разностное уравнение

Уравнения (1.13) недостаточно для определения двух неизвестных величин, поэтому поступим следующим образом [2]. Представим значение температуры в узле сетки N через разложение в ряд Тейлора в окрестности точки N -1:

т = Г»,-і +Аг, ||Т

(Агі)

2 ґд 2Г Л

2

. дг ,

\иГ У»,-1

Но из уравнения (1.1)

д2Т 1 дТ 1 дТ

тогда

дг2

дТ

дт г дг

ты, - ты, -і +*41 дг

(Аг1 )2 (1 дТ 1 дТ

2

дт г дг

Запишем полученное уравнение в конечных разностях:

грП + 1 грП

,-1~ ,-1

2а

Ат

(Агі)

(грП +1 'Т’П +1 \

N-11

а 2а ^ грП + 1 грП + 1 ^ ~тм,-2

^ Гм,-1 АГ1 1 2Аг,

(1.14)

Из уравнений (1.13) и (1.14) получим

т„+1 _ АГ1 2г«, -1 т„+1 + ^ _ Аг + 2г ^-2 +

Ш1 ^ 2,Ы1 -1

К,-1 (А^1 )2 + 4а^,-1Ат „ + 1

гг>п+ і ____

ТН, +1 ~ "

аАт (А Гі + 2гм,-1)

2гн,-і (Агі У

аАт (Аг1 + 2гК]-і) М' 1 Я,

- (1.15)

гг<п

N і + 1 * л

1 а А г + Я_

■Т + ±

N. т А Л N1 + 2

1 а А г + Я '

тп::2. (1.16)

п+1

N1 +1

N-1

Г

к

V

а

+

а

N-1

N -1

+

+

= а

Для определения температуры в узле N + N воспользуемся граничным условием (1.7):

7'П

-\Т

(1.17)

Для определения шага по времени Аг исследуем уравнение на устойчивость методом гармоник [3]. Представим решение разностной задачи в узле сетки в виде Т” = Япвл<р и подставим в уравнение (1.10). Получим, что схема устойчива при

V а г / ( Аг2 /'

Аг < тіп

2а„

(1.18)

Для решения поставленной задачи была написана программа на языке С# для расчета распределения температуры в стенке массива и воздухе.

В результате для ряда входных данных (Я1, Я2, и, То1, Т02, Дг, срг, рг, Яс, сс, рс) мы получили следующее графическое решение:

пренебрегаем, тогда, поскольку Т = Т(г,т) , уравнение, описывающее изменение температуры в газе, примет вид

дТ +д ( т) — д2т + Qv

-----I----(иТ) — а —— +----------.

дт дz дz ср

(2.1)

Для определения объемной плотности внутренних источников теплоты <2у выделим участок канала

длиной I, пусть й - диаметр канала, тогда тепловой поток от стенки канала к воздуху будет

Ч = « (Тст - т)п^.

Эквивалентный тепловой поток от объемных источников тепла

п ,

Я = -

Из этих равенств получаем, что в одномерной постановке

От =

2а (Тст - Т)

Я,

Дополним уравнение (2.1) уравнением неразрывности

^ + — (ир) = 0, (2.2)

дт ^ ’

уравнением движения

ди ди 1 др

— + и — +-------- = 0 (2.3)

дт дг р дг

и уравнением состояния

р = рЯТ . (2.4)

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

2. Решение задачи нестационарного процесса сопряженного теплообмена между рудничным воздухом и горным массивом в двумерной постановке

Пусть по цилиндрическому каналу радиуса & и длиной I движется однородный однонаправленный поток газа (воздуха) со скоростью и = иг (z,т). Толщина стенки канала &2 _ & , причем >> &

(рис. 1). Температура газа в начальный момент времени т = 0 во всем канале равна Ті, температура стенки в начальный момент времени составляет Т2, причем Т2фТі. Температура на внешней поверхности стенки канала в любой момент времени т > 0 остается постоянной и равной Т2.

Так как расчетная область обладает осевой симметрией, задачу будем решать в цилиндрических координатах. Ось Oz цилиндрической системы координат направим по оси канала (рис. 1). Поперечными течениями изменения температуры в газе

Здесь Я = 287

Дж

газовая постоянная воздуха,

кг ■ К

Т - абсолютная температура (в К).

Изменение нестационарного поля температур в стенке Тс = Тс(г^,т) описывается уравнением

дТс

дт ~ а°

-------2—+-----------+----2~

дг г дг дz

(2.5)

Приведем систему уравнений (2.1) —(2.5) к безразмерному виду. В качестве обезразмеривающих па*

раметров возьмем характерную плотность р , характерное давление р*, характерную длину Г , характерную температуру Т*, характерное время * / * р _ т = —, характерную скорость и = . — . При таи ]/ р

ком выборе обезразмеривающих параметров уравнения (2.2) и (2.3) сохраняют свой вид. Уравнение

(2.1) примет вид

д? д ( т\ - - д2т ^ -----+---(иГ) — а —— Л---------,

дт дz дz

Р

уравнение (2.5) примет вид

дТс _ дт ~ °с

д 2Тс 1 дТс д 2Тс дг2 г дг дz2

\

(2.6)

(2.7)

а уравнение (2.4) примет вид

V

р = рЯТ

(2.8)

* *

Здесь Я =—— Я — безразмерная газовая постоян-

ная, аг = —тї - безразмерная температуропровод-и 1

ность газа, ас = -**— безразмерная температуро-и 1

проводность

стенки,

От =

1* От

* * * Тирс

безразмерная объемная плотность внутренних источников теплоты.

При расчетах принималось: р* = 101325Па ,

*

Т* = 273К , /* = 2Я1, р* = -Аг .

1 ЯГ*

Систему уравнений (2.2), (2.3), (2.4), (2.5) дополним начальными условиями:

Т(2,0)|= Т, (г) ,

и)|о*= и, (г) ,

Р ( г’° )| о* * = Р, ( г) >

Тс (г, г, 0 )| ^ г ^ ^ = Тс1 (г,1) и граничными условиями:

Тс(г,2,т)|г=К21о<г<1 = Тс2 (г,т) ,

Тс(г.О.т)|К1 <^ = Тс3 (г,т) ,

Тс(г,1,т)|^^ = Тс4 (г,т),

а

(Т - Тс)

=; дТС|

Г-й7,0<2-І ^ дг г =К!,°-г~1 '

и(0,т) = и2 (т), и(1,т) = и3 (т),

Р (°-Т) = Р2 (Т) , Р (/,Т) = Рз (Т) Т(0,т) = Т2 (т), Т(1,т) = Т3 (т).

У

Введем вектор X = <

тогда систему уравнений

(2.2), (2.3), (2.6), описывающих изменение параметров воздуха во времени и по длине выработки, можно записать в виде

дХ дХ _д2 X -

---+ А---+ В—— = С ,

дт дz дz

где матрицы А =

V

Ґ \

0

0

ё

ґ и Р 0 ^ ( 0 0 0 '

ЯТ

и Я , в = 0 0 0

Р 0 0 ~—р

, 0 Т и V

а вектор С =

Р

Разностную схему, позволяющую получать по известным параметрам на п-м временном слое параметры воздуха в узлах сетки на (п+1)-м временном слое, запишем в виде

х;+1 = х; - л; ^ (х” - х;_1) -

-в; Ат (X;+1 - 2Х; - X; )+с; ат

или, переходя от векторной формні записи к компонентной,

Ат ^ п п \ п А.т і

—( р, - р,-1)- р,

п + 1 п п ( П п \ п Ат ( п п \

рі - рі - и, Т7(р, - р,-1)- р, Т7(и, - и,-1),

ЯТп Ат / \

п + 1 п ] 1 / п ”4

и , = и ,-----;-----(р; - р)_і ) -

)

Р” Аz

п Ат / ” п \ гл Ат (грп грп \

-иу ^7()-Я~А<Т'"Т'-і )•

т.„ = г; - т, - Пі).

р

+а ^ -2т?+т^.-1)

Здесь Ат — шаг по времени, Az — шаг по про-

^ ТП гп+1

странственной координате 2, и — значения

сеточной функции / в 7-м узле на п-м и п+1-м временных слоях.

Конечно-разностный аналог уравнения (2.7) запишем в виде

7СП7 = ТсП , + Ас Ат

'к,

' ТсП+1, у - 2ТСкИ ; + Тс1_, у Аг2

7 тсп - Тсп Тси - 2ТСп + Тсп + 1 тск+и тск-1,] + тск,] +1 2тск,ут ТСк,]-1

2Аг

Аг

Здесь Аг — шаг по пространственной координате г, Тспк ] и ТсП+; — значения температуры стенки в к-, ]-

узле сетки на п-м и п+1 -м временных слоях, гк = Я1 + k - А г .

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

Ат = тіп

А2-Аг2

тт — і и

П Аг - 2аг 2ас (Аг2 +Аг2)

3. Архитектура суперкомпьютера «ПГУ-Тесла»

Система «ПГУ-Тесла» — высокопроизводительный многопроцессорный вычислительный комплекс с гибридной архитектурой, произведенный ведущей российской суперкомпьютерной компанией «Т-Платформы».

к

Ґ

\

«ПГУ-Тесла» позволяет параллельно решать до 24Q задач на вычислительных ядрах процессоров Intel Xeon 567Q и дополнительно до 5376 подзадач на вычислительных ядрах процессоров Nvidia Tesla S2Q5Q.

Высокоскоростная сеть QDR Infiniband дает возможность в десятки раз быстрее обмениваться информацией между задачами суперкомпьютера, чем обычная сеть Gigabit Ethernet. Производительность системы:

9,0 Терафлопс - пиковая,

4,9 Терафлопс - на тесте Linpack (большие задачи линейной алгебры).

Суммарный объем оперативной памяти -96Q гигабайт.

4. Эффективность применения распараллеливания на GPU

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

В таблице приведены результаты решения задачи на стационарном ноут6укє(Ш:є1 Dual-Core T42QQ), процессоре Intel Xeon 567Q (ПГУ-Тесла), на видеокарте Nvidia Tesla S2Q5Q (непосредственный алгоритм на CUDA C, а также последовательный код с использованием директивы OpenACC). Задача решалась на сетке 1Q24 х 1025 в области [Q, 1QQ] х [Q, 20] , было выполнено 50000 шагов по времени.

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

Устройство Время решения Ускорение

Intel Dual-Core T4200 34 ч. 28 мин. 11 сек. -

Intel Xeon 5670 21 ч. 09 мин. 44 сек. l.65x

Nvidia Tesla S2050 (OpenACC) 03 ч. 37 мин. 48 сек. 9.4x

Nvidia Tesla S2050 (CUDA C) 00 ч. 46 мин. 35 сек. 45х

(ОрепАСС) 52050(си0д

ч

Сравнение времени решения задачи на СРи, GPU

5. Результаты расчетов

Проведены многочисленные расчеты с различными параметрами.

Для примера приведем результаты расчета программы на сетке 1024 х 1025 в области [0, 100] х [0, 8].

Прошло времени:

56мин. 2 ч. 57 мин.

8 м.

Т

Lz

б ч. 51 мин. 13 ч. 37мин.

I г I—2

6. Заключение

Работа направлена на численное моделирование нестационарного процесса сопряженного теплообмена между горным массивом и рудничным воздухом с использованием высокопроизводительной вычислительной системы «ПГУ-Тесла».

В процессе решения поставленной задачи произведена оценка предыдущих работ по данной тема-

тике [1—3]. Интересен тот факт, что практически все их авторы [1—3] основывают свои расчеты на одной и той же приближенной физико-

математической модели теплообмена, разработанной А.Н. Щербанём [1]. Однако область применения данной физико-математической модели теплообмена ограничена размерами временных

интервалов. Хорошая точность обеспечивается при расчете теплообменных процессов продолжительностью от 2 сут. и более. Для описания быстрых процессов, например, внезапного возгорания длительностью в часы, минуты, данная модель неприменима по причине больших погрешностей.

В нашей работе предложен другой подход к построению физико-математической модели теплообмена между горным массивом и рудничным воздухом.

В процессе решения поставленной задачи была реализована программа для моделирования тепло-переноса в горном массиве и рудничном воздухе с

применением технологии NVIDIA CUDA, а также с использованием директивы OpenACC.

Результаты работы будут использованы в программном комплексе «Аэросеть» Г орного института УрО РАН.

Список литературы

1. Исаченко В. П., Осипова В. А., Сукомел А. С. Теплопередача. М.: Энергия, 1975. С. 142—149.

2. Луканин В.Н., Шатров М.Г., Камфер Г.М. Теплотехника: учебник для вузов / под ред. В.Н. Луканина. М.: Высш. шк., 2000. С.229—241.

3. Михеев М. А. Основы теплопередачи. М.: Гос-энергоиздат, 1956. С.92—126.

Simulation of non-stationary process of conjugate heat mountain and mine air using technology computing on GPUs NVIDIA CUDA

A. R. Kutcev

Perm State National Research University, Bukirev St., 15, 614990 Perm

Existing models of heat transfer between the mine air and mountains YaV-lyayutsya overly simplistic and does not take into account a number of important physical processes in mines. Their main drawback is the lack of conjugation in the calculation of temperature fields and the array of air, when the temperature field in the array is unchanged. In this work, the conjugate heat transfer is modeled by solving the more complete equations of heat transfer.

The solution to this problem will allow for prediction of the thermal regime of deep ore defenders and to develop effective measures to normalize the climate in the work areas of underground mining.

As the computational domain is taken excavation cylindrical shape, with the problem is axisymmet-ric. For each medium is solved numerically unsteady heat equation, are coordinated by the Newton boundary condition for the heat flux. The results of the calculation are presented in the form of diagrams of the air temperature and the array, depending on the time. We study the system's output on a stationary thermal conditions.

Keywords: Transient process, NVIDIA CUDA, conjugate heat transfer, mining, mine air.