Numerical analysis of heat gain to geothermal heat exchangers with moisture evaporation in the structure

Full Text

Введение

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

Например, Исландия [4, 5], где отопление домохозяйств за счет геотермальных ресурсов приближается к 100 %, является в рассматриваемой области признанным лидером. Российская Федерация также обладает достаточными запасами геотермальной энергии и по оценкам, приведенным в [6], они в 8–12 раз превышают потенциал всех углеводородных ресурсов.

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

Основным недостатком аналитического описания процессов переноса в рассматриваемых системах является наличие достаточно большого количества допущений, без введения некоторых из них получить аналитическое решение невозможно. Можно выделить следующие главные ограничивающие допущения, при которых решаются задачи: однонаправленность теплового потока [7], использование граничных условий первого рода [8] или линейного источника тепла [9] на границах области расчета, возможность сегментирования конечного линейного источника тепла [10] и его стратификации [11].

Численное моделирование базируется чаще всего на использовании коммерческих программных комплексов: Feflow [12, 13], OpenGeoSys [14], Ansys Fluent [15, 16] и Comsol Multiphysics [17], в основу которых заложен метод конечных элементов. Основная трудность в случае использования такого программного обеспечения состоит в необходимости использования большого количества элементов пространственной сетки, а это обстоятельство неизбежно приводит к увеличению длительности вычислений. В исследовании [18] показано, что время расчета одного отопительного периода при эксплуатации геотермальной скважины на персональном компьютере с процессором 2,2 ГГц и 64 ГБ оперативной памяти составляет 143 часа.

Использование оригинальных программных кодов и, в частности, метода конечных разностей позволяет существенно сократить продолжительность расчетов до адекватных значений (несколько минут). В [19, 20] показано, что скорость расчетов рассматриваемых систем можно увеличить на порядки по сравнению с коммерческим программным обеспечением.

Отдельной областью в сфере моделирования геотермальных технологий [1–3] является исследование тепловых режимов геотермальных теплообменников (ГТ). Конструктивно наиболее типичный ГТ представляет собой теплообменник типа «труба в трубе», установленный в скважине (рис. 1). Необходимость обеспечения надежного теплового контакта между ГТ и окружающей его средой приводит к использованию разнообразных засыпок. В качестве засыпок используются различные материалы: от обычного цемента [21] до материалов с фазовыми переходами [22]. Однако использование этих материалов заметно повышает и без того существенные капитальные затраты при строительстве скважин с ГТ. В [23] было предложено использовать в качестве засыпки распространенный и недорогой материал – увлажненный песок. Использование увлажненного песка в качестве засыпки для ГТ неизбежно будет связано с наличием фазовых переходов и изменением механизмов теплопереноса в конструкции геотермального теплообменника.

 

Рис. 1. Схема скважины с ГТ: 1 – внутренняя труба; 2 – кольцевой канал; 3 – засыпка; 4 – обсадная колонна; 5 – горячие породы

Fig. 1. Diagram of a well with a geothermal heat exchanger: 1 – inner pipe; 2 – ring channel; 3 – backfill; 4 – casing; 5 – hot rock

 

Анализ научной литературы по исследуемой проблематике [1–23] показал, что исследование теплопереноса в конструкциях и зонах размещения ГТ с учетом взаимосвязи характеристик засыпок, режимных параметров и фазовых переходов до настоящего времени не выполнено. Целью данной работы является исследование влияния интенсификации теплоподвода к ГТ за счет изменения теплопроводности в результате увлажнения засыпки в элементах его конструкции и наличия фазовых переходов.

Постановка задачи

Прототипом рассматриваемой конструкции ГТ является реальный объект, геометрические и физические параметры которого описаны в [22]. Обсадная колонна, выполняющая в [22] роль несущей конструкции, выполнена из высокопрочного бетона, а кольцевой канал ГТ – из стали. На рис. 2 приведена схема области решения рассматриваемой задачи.

 

Рис. 2. Схема области решения: I – внешняя труба кольцевого канала ГТ; II – засыпка; III – обсадная колонна; IV – горячие породы

Fig. 2. Diagram of the solution area: I – outer ring channel pipe; II – backfill; III – casing; IV – hot rock

 

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

Влага, наполняющая поры песчаной засыпки II, на границе контакта (R3) с разогретой конструкцией обсадной колонны III испаряется. Образовавшийся пар диффундирует в зону меньших концентраций. На поверхности R2 происходит конденсация пара с выделением тепла (рис. 2).

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

Теплофизические характеристики материалов, веществ и их компонентов являются постоянными и известными величинами.

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

В местах соприкосновения слоев (рис. 2) выполняются условия идеального теплового контакта.

Не рассматривается перенос тепла в энергоносителе, циркулирующем в ГТ, а также процессы массопереноса в засыпке и горячей породе (рис. 2).

Считается, что вдали от скважины поддерживается неизменная температура пород.

Математическая модель

Уравнения теплопроводности и соответствующие им краевые условия для рассматриваемой области решения (рис. 2) имеют следующий вид:

$$\begin{gathered} \tau \ge \mathrm{0,}{R}_1\le r\le R_5; \\ {c}_i{\text{ρ}}_i\frac{\partial T_i}{\partial \tau }={\text{λ}}_i\left(\frac{{\partial }^2{T}_i}{\partial r^2}+\frac{1}{r}\frac{\partial T_i}{\partial r}\right);i=\text{I}-\text{IV;} \end{gathered}$$ (1)

$\tau =\mathrm{0,}{R}_1\le r\le R_5;T_i=T_0=\text{const};i=\text{I}-\text{IV;}$ (2)

$\tau >\mathrm{0,}r=R_1;\alpha \left(T_{r=R_1}-T_{\text{in}}\right)=-{\text{λ}}_{\text{I}}\frac{\partial T_{\text{I}}}{\partial r};$ (3)

$\tau >\mathrm{0,}r=R_2;{\text{\hspace{0.33em}λ}}_{\text{I}}\frac{\partial T_1}{\partial r}={{\text{λ}}_{\text{I}}}_{\text{I}}\frac{\partial T_2}{\partial r}+Jq;T_1=T_2;$ (4)

$\tau >\mathrm{0,}r=R_3;{{\text{\hspace{0.33em}λ}}_{\text{I}}}_{\text{I}}\frac{\partial T_2}{\partial r}={\text{λ}}_{\text{I}}{{}_{\text{I}}}_{\text{I}}\frac{\partial T_3}{\partial r};T_2=T_3;$ (5)

$\tau >\mathrm{0,}r=R_4;{\text{\hspace{0.33em}λ}}_{\text{I}}{{}_{\text{I}}}_{\text{I}}\frac{\partial T_3}{\partial r}={\text{λ}}_{\text{IV}}\frac{\partial T_4}{\partial r};T_3=T_4;$ (6)

$\tau >\mathrm{0,}r=R_5;T_{\text{5}}=T_{ex}=\text{const.}$ (7)

Теплофизические свойства песчаной засыпки при известных плотности ρ и объемной влажности W вычислялись из следующих соотношений [23]:

${\lambda }_{\text{II}}=-\mathrm{1,337}+\mathrm{0,00125}\text{ρ}+\mathrm{0,01}W\text{;}$ (8)

$c_{\text{II}}=-\mathrm{0,018}+\mathrm{0,0009}\text{ρ}+\mathrm{0,031}W.$ (9)

Для расчета коэффициентов теплоотдачи в канале кольцевого поперечного сечения использовалось широко известное уравнение [24]:

$$\begin{gathered} \text{Nu}=\mathrm{0,017}{\mathrm{Re}}^{\mathrm{0,8}}{\mathrm{Pr}}_{\text{e}}^{\mathrm{0,4}}{\left(\frac{{\mathrm{Pr}}_{\text{e}}}{{\mathrm{Pr}}_{\text{w}}}\right)}^{\mathrm{0,25}}{\left(\frac{d_2}{d_1}\right)}^{\mathrm{0,18}}, \\ \text{Nu}=\frac{\text{α}d}{{\text{λ}}_{\text{e}}};\mathrm{Re}=\frac{Vd}{{\text{ν}}_{\text{e}}};d=d_2-d_1. \end{gathered}$$

Коэффициент аккомодации, скорость и теплота фазового перехода вычислялись из соотношений [25]:

$$\begin{gathered} a=\frac{\mathrm{0,059}}{P_{\text{s}}^{{}^{\mathrm{0,5}}}}; \\ J=\frac{a\left(P_{\text{s}}-P_{\text{p}}\right)}{{\left(\frac{2\pi R_g}{M}{\left.T\right|}_{r=R_2}\right)}^{0.5}}; \\ q=\mathrm{2500,64}-\mathrm{2,369}{\left.T\right|}_{r=R_2}. \end{gathered}$$

Теплопритоки Q без учета термического сопротивления металлической стенки трубы ГТ находились из выражения:

$Q=-{\text{λ}}_{\text{II}}\frac{\partial {\left.T_2\right|}_{r=R_2}}{\partial r}\frac{F}{l}+\frac{JqF}{l}.$

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

Обозначения: λ – коэффициент теплопроводности, Вт/(м·К); α – коэффициент теплоотдачи, Вт/(м²·К); ρ – плотность, кг/м³; ν – кинематическая вязкость, м²/с; τ – время, с; R – граница области расчета, м; T – температура, К; c – теплоемкость, Дж/(кг·К); r – координата, м; V – скорость, м²/с; W – объемная влажность, %, d – эквивалентный диаметр, м; d2, d1 – внешний и внутренний диаметры кольцевого канала, м; Nu, Re, Pr – числа Нуссельта, Рейнольдса и Прандтля; π=3,14; a – коэффициент аккомодации; J – скорость испарения, кг/(м²·с); М – молекулярная масса паров, кг/моль; P – давление, Па; q – теплота фазового перехода, Дж/кг; Rg – газовая постоянная, Дж/(кмоль·К), Q – теплопритоки, Вт/м; l – единичная длинна, м; F=2πR2l – площадь поверхности, м².

Индексы: 0 – начальный момент времени; 1, 2, 3, 4, 5 – номера границ областей расчета (рис. 2); in – внутренний; ex – наружный; e – жидкость (энергоноситель); w – стенка; p – парциальный; s – насыщение; I, II, III, IV – номера областей расчета (рис. 2).

Метод решения и исходные данные

Решение задачи (1)–(7) получено методом конечных разностей. Использовалась неявная разностная схема и алгоритм прогонки. Шаг по координате составлял от 1 до 10 мм. Геометрические параметры (R1=0,05 м; R2=0,055 м; R3=0,14 м; R4=0,25 м) соответствовали конструкции скважины, описанной в [22], а R5=10 м. Внешний и внутренний диаметры кольцевого канала составляли: d1=2·R1; d2=0,06 м.

Начальная температура Т0 в рассматриваемой области решения принималась равной температуре разогретых пород Тex=373,15 К. Температура прокачиваемого энергоносителя составляла Тin=278,15 К. Объемная влажность песчаной засыпки W варьировалась от 5 до 25 % и была ограничена открытой пористостью. Расход энергоносителя задавался от 0,004 до 0,04 м³/с, что соответствует скоростям движения по кольцевому каналу V=0,1–1 м/с.

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

 

Таблица 1. Теплофизические характеристики

Table 1. Thermophysical characteristics

Характеристика Characteristic	λ, Вт/(м·К) W/(m·K)	c, Дж/(кг·К) J/(kg·K)	ρ, кг/м3 kg/m3
Кольцевой канал Ring channel	57,5	466	7860
Засыпка/Backfill	Расчет по (8) Calculation by (8)	Расчет по (9) Calculation by (9)	1900
Цемент/Cement [22]	1,78	800	2490
Горячие породы Hot rocks [22]	1,3	775	1990

 

Результаты численного моделирования

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

 

Таблица 2. Изменение теплопритоков к ГТ с учетом испарения при увеличении W от 5 до 25 %, %

Table 2. Change of heat gain to the geothermal heat exchanger taking into account evaporation with increasing W from 5 to 25%, %

τ, мес. month	V, м/с/m/s
	0,10	0,25	0,50	1,00
0,03	7,12	7,21	7,49	7,79
0,1	5,58	5,64	5,67	5,69
0,5	4,32	4,36	4,38	4,39
1	3,93	3,97	4,00	4,09
2	3,59	3,60	3,61	3,64
3	3,42	3,44	3,46	3,49
4	3,31	3,33	3,35	3,37
5	3,23	3,25	3,26	3,27
6	3,16	3,18	3,19	3,20

 

Таблица 3. Изменение теплопритоков к ГТ с учетом испарения при увеличении V от 0,1 до 1,0 м/с

Table 3. Change of heat gain to the geothermal heat exchanger taking into account evaporation with increasing V from 0,1 to 1,0 m/s

τ, мес. month	W, %
	25	15	5
0,03	2,47	2,40	2,31
0,1	1,65	1,61	1,57
0,5	1,47	1,43	1,40
1	1,34	1,31	1,27
2	1,22	1,19	1,17
3	1,05	1,03	1,01
4	1,12	1,11	1,08
5	1,14	1,13	1,10
6	0,96	0,94	0,93

 

Таблица 4. Рост теплопритоков к ГТ за счет испарения, %

Table 4. Increase of heat gain to the geothermal heat exchanger due to evaporation, %

W, %	V, м/с/m/s
	0,10	0,25	0,50	1,00
5	22,02	21,89	21,91	21,84
15	21,69	21,57	21,58	21,52
25	21,41	21,29	21,30	21,24

 

Рис. 3. Снижение теплопритоков к ГТ: 1 – W=25 %; 2 – W=15 %; 3 – W=5 %

Fig. 3. Reducing heat gain to geothermal heat exchanger: 1 – W=25%; 2 – W=15%; 3 – W=5%

 

Рис. 3 в качестве примера иллюстрирует типичную картину снижения теплопритоков Q к ГТ во время эксплуатации при фиксированной скорости движения энергоносителя по кольцевому каналу (V=1,0 м/с), различных значениях объемной влажности песчаной засыпки W и учете испарения влаги в засыпке.

Изменение величин теплопритоков к ГТ свидетельствует об их ожидаемом снижении с увеличением длительности эксплуатации из-за охлаждения горячих пород в непосредственной близости от скважины (рис. 1) и закономерном увеличении при росте W, сопровождающемся изменением теплофизических характеристик засыпки в соответствии с формулами (8) и (9). Время эксплуатации ГТ (6 месяцев) выбрано исходя из соответствия типичной для РФ продолжительности отопительного периода.

Анализ нестационарности процессов теплопереноса в рассматриваемой системе в период установлении стационарного температурного поля в горячих породах свидетельствует о её существенном влиянии на уровень теплопритоков к ГТ. За рассматриваемый период эксплуатации (6 месяцев) теплопритоки к ГТ снижаются практически в 2 раза по сравнению с первоначальным уровнем. Здесь следует отметить, что результаты моделирования указывают на резкое снижение тепловых потоков в первые дни работы ГТ (аналогичные выводы сделаны в [22] при исследовании материалов с фазовыми переходами). Это объясняется интенсивным охлаждением металлического корпуса ГТ (рис. 1) из-за его сравнительно высокой теплопроводности (табл. 1). Этот результат имеет конкретное практическое применение, поскольку позволяет обоснованно выбирать пути регулировки расхода энергоносителя для выравнивания теплосъема от ГТ во времени.

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

Исследование влияния влажности песчаной засыпки (W=5–25 %) на теплопритоки в конструкции ГТ позволило сделать вывод об увеличении теплопритоков до 7,8 % в начальный период эксплуатации ГТ и до 3,2 % к концу этого периода (табл. 2). Более высокий рост теплопритоков в начальный промежуток времени объясняется теми же обстоятельствами, что и при анализе нестационарности процессов переноса в рассматриваемой системе. Несмотря на достаточно скромное увеличение тепловых потоков в системе геотермальной скважины с ГТ из-за изменения объемной влажности засыпки, суммарный рост теплопритоков в рассматриваемой системе при глубине скважины 25 м за весь период эксплуатации (6 месяцев) может приводить к дополнительной аккумуляции в 15 Гкал и более.

Анализ влияние интенсивности теплоотдачи в кольцевом канале на изменение теплопритоков к ГТ (табл. 3) позволяет сделать вывод о не значительном росте тепловых потоков в рассматриваемой системе (около 1–2,5 %). Это объясняется достаточно быстрым установлением постоянной температуры на внешней трубе кольцевого канала ГТ (рис. 2). Следовательно, при решении задач, подобных задаче (1)–(7), для исследования тепловых режимов ГТ можно обоснованно вместо выражения (3) использовать граничные условия первого рода.

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

В табл. обобщены результаты численного анализа влияния наличия испарения влаги в конструкции ГТ на рост теплопритоков к нему в зависимости от W и V в конце срока эксплуатации (6 мес.). Анализ содержания табл. 4 свидетельствует о существенном вкладе процесса испарения влаги в засыпке ГТ на интенсификацию теплопритоков к нему. Рост теплопритоков к ГТ с учетом наличия испарения составляет около 22 % по сравнению с аналогичным исследованием без учета этого эффекта.

Адекватность результатов численного моделирования следует из проверок используемых методов решения задачи (1)–(7) на сходимость и устойчивость, а также подтверждается численным сопоставлением с известными данными о работе геотермальных скважин с ГТ [21, 22]. В [21, 22] линейные плотности теплового потока составляют 60–300 Вт/м, что хорошо согласуется с результатами данной работы (рис. 3).

Заключение

Установлено, что увеличение влажности песчаной засыпки ГТ приводит к росту теплопритоков на 3,2–7,8 %.

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

Показано, что вклад эффекта испарения в теплопритоки к геотермальному теплообменнику составляет около 22 %.

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

About the authors

Viacheslav Yu. Polovnikov

National Research Tomsk Polytechnic University

Author for correspondence.
Email: polovnikov@tpu.ru
ORCID iD: 0000-0001-5504-0411

Dr. Sc., Professor

Russian Federation, 30, Lenin avenue, Tomsk, 634050

Sofia D. Shelemekhova

National Research Tomsk Polytechnic University

Email: shelemehovaaa@mail.ru
ORCID iD: 0000-0002-9898-8867

Postgraduate Student

Russian Federation, 30, Lenin avenue, Tomsk, 634050

Egor V. Lyubivy

National Research Tomsk Polytechnic University

Email: evl27@tpu.ru
ORCID iD: 0009-0009-8908-2041

Postgraduate Student

Russian Federation, 30, Lenin avenue, Tomsk, 634050

References

Supplementary files