Тепломассообмен в приповерхностных геотермальных системах Текст научной статьи по специальности «Строительство и архитектура»

УДК 621.311.1

Э.И. БОГУСЛАВСКИЙ, д-р техн. наук, профессор, boguslEI@yandex.ru Н.Н. СМИРНОВА, канд. физ-мат. наук, доцент, ninasmirnat@yandex.ru С.В. ЕГОРОВ, канд. физ-мат. наук, доцент, s.v.egorov@mail.ru

Санкт-Петербургский государственный горный институт (технический университет) E.I. BOGUSLAVSKIY, Dr. eng. sc.,professor, boguslEI@yandex.ru

N.N. SMIRNOVA, PhD in phys. and math. sc., associate professor, ninasmirnat@yandex. ru S.V. EGOROV, PhD in р.hys. and math., associate professor, s.v.egorov@mail.ru Saint Petersburg State Mining Institute (Technical University)

ТЕПЛОМАССООБМЕН В ПРИПОВЕРХНОСТНЫХ ГЕОТЕРМАЛЬНЫХ СИСТЕМАХ

Приведены результаты математического моделирования кондуктивного тепломассообмена в приповерхностных геотермальных системах. Теплообменник системы описывается как вертикальный канал в непроницаемом породном массиве. Для описания температуры пород используется уравнение теплопроводности в цилиндрических координатах. На границе пород и теплообменника выполняются условия III рода.

Ключевые слова: тепломассообмен, приповерхностные геотермальные системы, непроницаемый породный массив, численное моделирование.

HEATMASSTRANSFER TO NEARSURFACE GEOTHERMAL SYSTEMS

The paper is concerned with mathematical modeling of heat and mass transfer in subsurface geothermal systems. Heat exchanger is driven as a vertical channel in an impermeable rock massive. The temperature of rock is driven in accord with a usual thermal conductance equation on a cylinder coordinates basis. The boundary conditions between the rock and heat exchanger are given by the equality of a thermal flow.

Key words: heatmasstransfer, nearsurface geothermal systems, impermeable rock massive, numerical modeling.

Одним из основных ресурсосберегающих и экологически безопасных источников энергии является геотермальная. Технологии, позволяющие извлекать энергию земных недр, давно вызывают обоснованный интерес. Приповерхностные геотермальные установки - это наиболее доступный тип геотермальных систем. Они позволяют утилизировать энергию пород, залегающих на глубине не более 200 м [1, 3].

Физическая модель задачи о кондук-тивном теплообмене. Кондуктивный характер теплообмена подразумевает только потоки тепла, обусловленные разностью

температур пород и устройств образующих систему.

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

2R,

Рис. 1. Поперечное (а) и продольное (б) сечение приповерхностной геотермальной установки

1 - скважина; 2 - нагнетательная труба; 3 - засыпка между трубой и скважиной (песок - вода); Rl - радиус нагнетательной трубы; Тп - температура пород; d1 - диаметр нагнетательной трубы; d2 - диаметр скважины; и - скорость течения

теплоносителя

б

a

1

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

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

Математическая модель процесса кондуктивного теплообмена. Задача теплообмена пород и теплоносителя текущего в цилиндрическом канале математически описана, например, в работах [2, 6]. Рассмотрим цилиндрический канал, помещенный в среду с постоянной температурой Т = const. На внутренней границе с радиусом R, происходит теплообмен с непроницаемыми горными породами.

Теплообмен на границе цилиндра и среды подчиняется закону Ньютона (граничные условия III рода) [4, с.181].

Температура теплоносителя t(x, т) зависит от координаты x, отсчитываемой вдоль канала и от времени т. На входе в канал при х = 0 температура задается постоянной и равной t (0, т) = t0 = const. В начальный момент времени температура теплоносителя в канале равна температуре окружающей среды t (x,0) = Тп.

Внутри канала осуществляется конвективный перенос тепла в жидкости, в массиве горных пород - кондуктивный перенос тепла. В любой момент времени т в массиве горных пород имеется нестационарное распределение температуры Т = Т (x, r, т)по радиусу r и по глубине х, которое подлежит определению. В канале - нестационарное распределение температуры теплоносителя t = t (о, т) по глубине х, которое также подлежит определению. Зависимость температуры в канале от расстояния до центра канала не учитывается.

В этом случае дифференциальное уравнение для температуры теплоносителя в трубе и для температуры в массиве горных пород соответственно

5t(О, т) ^ f dt(о, т)

"ГМ"

о

ст

дх

а

i (t (о, т) - Т (R,, х, т));

и

дТ(х,г, т) А,2

дт

Р2С2

c,Pi

( д2Т(х,г,т) +1 дТ(х,г,т)Л дг2 r дг

_ 25

Санкт-Петербург. 2010

.Xo

.ro Г .r

Рис.2. Сетка математической численной модели процесса кондуктивного теплообмена

при начальном условии £ (0, т) = £0, Т(х, г,0) = Т и граничном условии

, дТ

— Л 2-

дг

= «i(t — Т )1

г=Ri

Здесь и - скорость течения теплоносителя по трубе; а - отношение площади поверхности теплоносителя к его объему, для цилиндра а = 2/Яь Я1 - радиус трубы; а1 -коэффициент теплоотдачи границы раздела теплоносителя и породы, □ который определяется геометрическими характеристиками системы и тепловыми свойствами засыпки; р1 и с1 - плотность и теплоемкость теплоносителя; т - время; £ - текущая температура теплоносителя; Т - текущая температура в массиве горных пород; х - координата вдоль трубы; г - текущий радиус; Л,2, с2 и р2 - теплопроводность, удельная теплоемкость и плотность пород в массиве соответственно.

Так как теплообмен между теплоносителем и массивом начинается с момента прихода гидродинамической волны в данную точку с координатой х, перейдем к но*

вым переменным х, т = т - х/и. В новых координатах система дифференциальных уравнений и граничных и начальных условий запишется в виде

и (х,т*) _Т{Ях,х,т*)]; (1)

дх

ciPi

дТ (х, г, т*)

дТ*

Л Гд2Т(х,г,т*) 1 дТ(х, г,т*)Л

Р2С2

дг2

- + — г

дг

Граничное условие III рода примет вид

. дТ

— Л 2-

дг

= a1(t — Т )|

г=R1

г = Ri

начальное условие

t(o,т*) = t0,Т(х,г,0) = Т .

(3)

(4)

Дифференциальные уравнения (1) и (2) с граничным условием (3) и начальным условием (4) согласовано решались численным методом. Согласно сетке численной задачи (рис.2) координата х отсчитывается вдоль трубы в вертикальном направлении, координата г - в горизонтальном направлении; I - номер узла сетки в направлении х; ] - номер узла сетки в радиальном направлении. Начальная точка оси х соответствует поверхности земли и обозначается х0, конечная точка хтах - концу трубы. Длина трубы равна длине нетеплоизолированного участка трубы теплообменника, внутренний радиус трубы г0 соответствует стенке трубы. Шаги сетки вдоль оси х и вдоль оси г обозначаются Ах и Дг соответственно и задаются независимо. Процедура перехода к конечным разностям аналогична той, что описана в работе [5, с.221].

Уравнения решаемой системы, будучи записаны в виде конечно-разностных схем, имеют вид

и!Сх±1Х)_!(хл1 =

Дх

а

ciPi

(i[t (х, тк *) — Т (х, ^ тк )]>

Т(х, г, тк+i)—Т(х-, г, тк) _ л2

Ат*

Р2С2

Т (х, j, тк)—2Т (х, г, тк)+Т (х, г—, тк)

Аг2

+

+

iТ (х, гтк)—Т (х, г, тк)

Аг

Граничное условие

х

X

г

-А

Аг

Т (хг, ^ Т-) - Т (хг, гo, Т*к )

1!'(х-, Тк) - Т (х, гo, ч)].

= а1

Начальные условия запишем в виде

'(х0, т*к) = при любых к и Т(хг, -, т0) = Т

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

Из первого и второго уравнений данной системы получим соответственно

'(хк)=

аАх

с1р1и

х а1 ^(х, Тк *) - Т (х, ^ тк)]+'(х, Тк *); (5)

Т (х, -, тк+1)=ААТ-х р2С2

Т (х-, 0+1, тк) - 2Т (х-, -, тк)+Т (х-, г-, тк)

Аг2

+

+ -

1Т (х, тк) - Т (х, , т- )

Аг

+

+Т (х, г,, тк).

Граничное условие примет вид

Т (X, ^ Тк) = Т (х, гl, Тк) + + а1 х- '(X, Тк) - Т (X, ^ Тк)].

(6)

(7)

Эти формулы позволяют моделировать временную динамику температурного поля пород, сопровождающую движение теплоносителя в трубе, и изменение температуры самого теплоносителя по мере его движения. Процесс непрерывного перемещения теплоносителя дискретизируется. Считаем, что теплоноситель двигается вдоль оси х дискретно с шагом Ах. При этом температура теплоносителя меняется согласно формуле (5), а сопутствующее изменение температуры стенки трубы - формуле (7). По формуле (6) рассчитывается изменение температуры в массиве пород во времени. Шаг по времени Ат берется равным интервалу про-

движения фронта теплоносителя на величину Ах. При скорости теплоносителя и имеем Ат* = Ах / е . При начальной температуре массива пород Тп можно по формуле (5) рассчитать изменение температуры теплоносителя при первоначальном прохождении гидравлического фронта по трубе. Очевидно, что эффективное время т*к будет одинаково для всех точек хг (г = 0, 1, ... гтах), причем гтах определяется числом шагов по сетке и задается произвольно. Поскольку т*к соответствует моменту прихода фронта теплоносителя и началу процесса теплообмена на уровне хг, логично обозначить его т0 . Тогда температура в трубе

' (X, т0) =

аАх

с1р1и

а1 х

\(^ Т0-) - Т(х0 , ^ Т0-)]+ '(х0 , Т0-) •

* = аАх I (х2, Т0) — а 1 х с1р1и

< !(х1, Т0-) - Т(х1, г0, Т0*)]+ '(х1, Т0-) ;

'(XN , Т0) =

аАх

-1

с1р1и

х \(хм-1, Т0-) - Т(хм-1, Г0, Т0-)]+ '(хм-1, Т0-) .

Здесь N - число узлов на оси х. Температура Т(хг, г0, т0) берется равной Тп при любом г. Зная распределение по глубине температуры в трубе, можно рассчитать распределение по глубине температуры стенки, используя формулу (7):

Т (X0, Гo, Т0) = Т (X0, гl, Т0) + +а1 V- '(xo, Т0) - Т (xo, гo, Т0) ];

А.,

Т (xl, гo, Т0) = Т (xl, гl, Т0)+

(8)

+а1'(xl, т0) - Т (xl, гo, Т0)]

А

2

Санкт-Петербург. 2010

х

х

г

Т (^, гo, т0) = Т (^, т0)+

± а Дг [(хы, т0) _ т(хм, ^ т0)].

Л 2

Зная изменение температуры стенки трубы, рассчитаем изменение температуры массива по формуле (6):

Т(Xo,г), х1) =

Л 2АТ Р2С2

Т (х0, j, Т0) — 2Т (х0, г), т0) + Т (х0, ji, т0) Аг2

+

+ -

i Т(xo, г+l, т0)—Т(xo, г, т0)

Аг

+Т (xo, г, т0);

Т(х1, г} , т*) = _Л2_Ат х Р2С2

Т (х1 , г+i, т0)—2Т ( х1 , г, т0)+Т ( х1 , j i, т0) + (9)

Аг

+

i Т(xl, г+l, т0)—Т(xl, г, т0)

Аг

Т (%, г, т*) =

+ Т(xl, г, т0);

Л 2АТ Р2С2

Т(XN , г +1, т0) — 2Т(^ , Г} , т0) + Т(XN , j, т0)

+

Аг2

i Т(^, гj+l, т0)—Т(^, гj, т0)

+

Аг

+Т (^, г, т0).

Здесь значение индекса ] определяется числом шагов по оси г.

В итоге мы получаем распределение температуры в массиве пород на 1-м шаге по времени. Шаг по времени Дт* = Дх / е . Отметим, что формулы (9) могут быть применены к любому узлу сетки в массиве по-

род кроме первого и последнего узла. Первый узел сетки соответствует стенке. Ее температура рассчитывается с помощью формул (8). В последнем узле температура приравнивается к температуре пород. Эти вычисления дают распределение температуры теплоносителя в трубе и температурное

поле в породе через время Дт* = Дх /е после прохождения фронта теплоносителя. Следующий шаг аналогичен описанному выше, но при этом на входе его используется распределение температуры, полученное на выходе предыдущего. Этот процесс повторяется произвольно заданное число раз.

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

3t (х, г, т) дт

Л1 Г д 2t (х, г, т) 1 дt (х, г, т)

РС

Л

дг1

дг

дТ (х, г, т) Л 2

3т

( д 2Т (х, г, т) +1 дТ (х, г, т)

Р2С2

дг2

дг

и граничными условиями на границе нагнетательная труба - порода

— Л

дí (х,г,т)

дг

= -Л.

дТ (х, г, т)

г=R11

дг

г=R11

Процедура перехода к численной схеме здесь аналогична описанной.

Результаты математического моделирования. На рис.3 приведены результаты моделирования температурного поля на глубине 200 м для скорости прокачки теплоносителя воды 0,2 м/с, что соответствует расходу теплоносителя 5,71 м3/ч. Все расчеты выполнены для теплообменников длиной

х

X

г

X

г

г

г

X

X

г

Тп, °С

10

Ъ?^ - C^5"ч

д о°++ А 2 ч

° 4 ч

о+ д 8 ч

о+

А

О

24 ч

о+ + 48 ч

0,1 0,2 0,3 0,4 0,5

Рис.3. Динамика температурного поля на глубине 200 м

0

г, м

200 м. Радиус нагнетательной трубы Я1 = = 0,05 м. Коэффициент теплоотдачи на границе нагнетательной трубы а1 = 21 Вт/(м -°С). Температура теплоносителя на входе в теплообменник = 1 °С, температура пород Тп = 10 °С. Теплофизические свойства теплоносителя и пород следующие:

Свойства Вт/(м'°С) с, Дж/(кг'°С) Р, кг/м3

Вода 0,57 4190 1000

Порода 2,5 1300 2400

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

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

теплообменника выходит количество теплоносителя с массой Ат = рАх2nR1, где р -

плотность теплоносителя; Ах - шаг сетки в направлении движения теплоносителя; Я1 -радиус нагнетательной трубы. Тепло, полученное этим количеством теплоносителя от пород, Аq = Атс(' - ), где с - удельная теплоемкость; ' - температура теплоносителя на выходе из теплообменника.

Сложив порции тепла, приобретенные теплоносителем на каждом шаге, найдем полное количество тепла, полученное теплоносителем от пород (рис.6). Полное количество тепла, отобранное от пород за двое суток, 1 • 109 Дж. Оцененная по этим данным средняя мощность теплоотбора равна примерно 6 кВт.

С другой стороны, зная конечное распределение температуры в породах, можно рассчитать количество тепла, отданное породой. Расчет дает нам 1,3 • 10 Дж. Видно, что тепло приобретенное теплоносителем, и тепло, отданное породами, в нашей модели хорошо совпадают.

_ 29

Санкт-Петербург. 2010

t, °С 3

2,5 2

1,5 1

10

20

30

40 Время, ч

Рис.4. Температура теплоносителя на выходе из теплообменника

Тп, °С

10

о о о

°1

0,2

0,4

Рис.5. Восстановление температурного поля на глубине 200 м после работы теплообменника (охлаждения) в течение 2 сут (1) через 1, 8 и 24 ч (2, 3 и 4 соответственно)

Q, 108 Дж 10 8 6 4 2

10

20

30

40

50 Время, ч

Рис.6. Тепло Q, отобранное от пород. Скорость прокачки теплоносителя 0,2 м/c; коэффициент теплообмена а! = 21 Вт/(м2-°С)

9

8

7

6

5

0

0

r, м

0

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

ЛИТЕРАТУРА

1. Богуславский Э.И. Оценка и перспективы теплового освоения приповерхностных (до глубины 200-300 м) геотермальных ресурсов России / Э.И.Богуславский А.А.Смыслов, Б.Н.Хахаев // Проблемы промышленной теплотехники. Киев, 2007.

2. Дядькин Ю.Д. Геотермальная теплофизика / Ю.Д.Дядькин, С.Г.Гендлер, Н.Н.Смирнова. СПб, 1993.

3. Литвиненко В.С. Геотермальные ресурсы России, технология и экономика их освоения // В.С.Литвиненко, Э.И.Богуславский // Геофонд-2003: Международный семинар. М., 2003.

4. Лыков А.В. Теория теплопроводности. М., 1967.

5. Ортега Дж. Введение в численные методы решения дифференциальных уравнений / Дж.Ортега, У.Пул. М., 1986.

6. Смирнова Н.Н. Нестационарный теплообмен при фильтрации в гетерогенных средах / Ин-т теплофизики СО АН. Новосибирск, 1990.

REFERENCES

1. Boguslavskii E.I., Smyslov A.A., Hahaev B.N. Estimation and prospects of thermal development of nearsurface (up to 200-300 m depth) geothermal resources of Russia // Industrial heating devices problems. Kiev, 2007.

2. Dyad'kin Yu.D., Gendler S.G., Smirnova N.N. Geothermal thermophysics. Saint Petersberg, 1993.

3. LitvinenkoV.S., Boguslavskii E.I. Geothermal resources of Russia, technology and economics of development // Geo-fond-2003: International meeting. Moscow, 2003.

4. Lykov A. V. Heat transfer theory. Moscow, 1967.

5. Ortega J., Pul U. Introduction in numerical methods applied to differential equation. Moscow, 1986.

6. Smirnova N.N. Nonstationary heat transfer under filtration in geterogeneous media. Novosibirsk, 1990.