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

Математические заметки СВФУ Октябрь—декабрь, 2014. Том 21, № 4

УДК 519.63

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ ТЕМПЕРАТУРНОГО ПОЛЯ

МНОГОЛЕТНЕМЕРЗЛЫХ ГРУНТОВ ПРИ ВОЗДЕЙСТВИИ ТРУБОПРОВОДОВ М. П. Акимов, П. Е. Захаров, О. И. Матвеева

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

Ключевые слова: фазовый переход, многолетнемерзлые грунты, метод конечных элементов.

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

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

© 2014 Акимов М. П., Захаров П. Е., Матвеева О. И.

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

1. Геометрия

Рис. 1. Схематическое отображение области.

Для учета неоднородности грунта и трубопроводов строим геометрию из нескольких подобластей. Обозначим пустую внутреннюю область трубопровода через (рис. 1), а остальные подобласти с материалами будем обозначать через ^, где а = 1,... ,М — номер трубы, М — количество трубопроводов, к = 1,..., Nа — номер материала, Nа — количество материалов.

Каждую подобласть для слоя грунта обозначим через Пв, где в = 1,... ^ — номер слоя, N — количество слоев грунта. Тогда область всех слоев обозначим через

N

п=и пв.

в=1

Введем обозначение расчетной области

м

п = пд \ и .

а=1

Для границ области введем следующие обозначения: Га — верхняя граница грунта, Г^ — нижняя граница грунта, Гд — боковые границы грунта,

м

г. = и

а=1

— внутренние стенки трубопроводов.

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

Рассматривается нелинейная нестационарная задача теплопроводности [3-6], которая описывается уравнением с фазовым переходом

Р(С + ^/Ж -сИу(Л£га<1Т) = °> 0<*<*е. (!)

где р — плотность, с — теплоемкость, ^ — индикатор ледосодержания, ш — водо-насыщенность, Ь — удельная теплота плавления-кристаллизации, Т — температура, 4 — время, Л — теплопроводность. Уравнение (1) дополняем граничными условиями: на боковых границах считаем, что поток тепла отсутствует:

дТ

А—=0, X е Г,, (2)

дп

на верхней границе области задаем граничное условие третьего рода, которое описывает теплообмен с воздухом [7, 8]:

дТ

-А— = аа(Т - Та), х € Га, (3)

дп

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

Т = Т* х е Г^, (4)

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

дТ

-А — = аг{Т-%), х€Г*. (5)

дп

Здесь используются обозначения: п — внешняя нормаль к границе, аа(£) — суммарный коэффициент теплообмена, Та(£) — температура воздуха, Т^ — температура грунта на нижней границе, а и Т4(х, 4) — коэффициент конвективного теплообмена и температура воды в трубопроводах.

Для замыкания системы уравнений задаем начальное распределение температуры

Т(х, 0) = Т0(х), х е о. (6)

Индикатор ледосодержания ^ определяем как кусочно-линейную функцию [9], которая зависит от текущей температуры:

Г 1, Т < Т* - Тд,

<р(Т) = \ Т*Х 7 / 1 1

{ 0, Т > Т* + Тд,

где Т* — средняя температура фазового перехода, Тд — значение, определяющее интервал фазового перехода.

Плотность, теплоемкость и теплопроводность записываем через индикатор ледосодержания и значения в талой и в мерзлой зонах

р(х,Т)= <^(Т)р-(х) + (1 - <^(Т)), с(х,Т) = <^(Т)с_(х) + (1 - <^(Т))с+ (х), Л(х,Т) = <^(Т)Л_(х) + (1 - <^(Т))Л+ (х), где х е О, р_(х),р+ (х) — плотность мерзлой и талой зоны, с_(х), с+ (х) — теплоемкость мерзлой и талой зоны, Л_ (х), Л+(х) — теплопроводность мерзлой и талой зоны.

Считаем, что в подобласти трубопроводов водонасыщенность равна нулю и теплофизические параметры талой и мерзлой зоны одинаковы:

ш(х) = 0, р_(х) = р+(х), с_(х) = с+(х), Л_(х) = Л+(х), х е

Температурное состояние грунта описывается системой уравнений (1)-(6).

3. Аппроксимация

Для аппроксимации задачи (1)—(6) по времени используем метод конечных разностей [10,11], а для аппроксимации по пространству — метод конечных элементов [12,13]. Вводим равномерную сетку для времени с шагом т:

Ш = {t = tn | tn = пт, п = 0,1,..., N, Nt = te}. (7)

Для основного уравнения (1) используем двухслойную неявную схему с линеаризацией коэффициентов по предыдущему временному слою [14-16]:

/ dinn \ Tn+1 — Tn РП с" + шЬ)--div(A™gradTn+1) = 0

¿Г ) т у ° ' ' (8)

x е П, п = 0,1,..., N - 1.

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

а(Гп+1,у) = /(у), V е V, (9)

где а(Гп+1,у) — билинейная форма

а(Тп+1,у) = ! + ^ шЬ^Тп+1у + т(Хп gradTrl+1, gradг^)^ сЬс

о

+/т<+1г"+,у>"а°>

Г„ ъ

а /(у) — линейная форма

/(у) = ! + ^ + J та"+1Т™+1г> ¿в + ^ тагТ?+1у <1,з. (И)

О Га Г4

Граничные условия (2), (3) учитываются в уравнении (9) при использовании формулы Гаусса — Остроградского, т. е. при переходе от объемного к поверхностному интегралу.

4. Численная реализация

Геометрия и триангуляция областей строилась в программе СМ8И [17]. Программа СМЯИ является генератором 2Б и 3Б сеток конечных элементов со встроенным редактором геометрии. Генерация сетки выполнена методом триангуляции Делоне. Численное решение задачи реализована с помощью библиотеки для научных вычислений ЕЕшСЭ [18]. Библиотека ЕЕшС8 предназначена для решения систем дифференциальных уравнений методом конечных элементов. В качестве конечных элементов выбраны лагранжевы элементы первого порядка.

Система трубопроводов состоит из двух отопительных трубопроводов и одного трубопровода с холодной водой. Геометрия области представляет собой

Рис. 2. Расчетная область и схема расположения трубопроводов.

Рис. 3. Сетка: 2623 узлов, 5054 ячеек.

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

Численное решение сравниваем с эталонным решением, которое получено при достаточно мелких сетках по пространству и времени. Количество узлов пространственной сетки для эталонного решения равен 125600 (249737 ячеек), а шаг по времени т = 3 часа. Кроме разности численного и эталонного решения рассматривается погрешность

£ =

|Т П _ Тп II

_1_е I I 1/2

\\1Wl, '

где ТП(х) — эталонное решение п-го временного слоя, скалярное произведение и норма определены следующим образом:

М) = У V2 ¿х, || V|| = (^)1/2.

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

Обозначение, ИГЭ-1 ИГЭ-2 ИГЭ-3 ИГЭ-4 ИГЭ-5

единица измер. (мерз./тал.) (мерз./тал.) (мерз./тал.) (мерз./тал.) (мерз./тал.)

а;, % 0,24 0,18 0,24 0,28 0,24

Г*, °С -0,64 -0,7 -1,06 -0,1 -0, 12

р, кг/м3 1710,0 1850,0 1760,0 1850,0 1900,0

с, Дж/(кг-°С) 1161,4/795,5 1107,2/811,0 1181,6/823,4 1136,4/827,9 1106,5/827,9

Л, Вт/(м °С) 1,80/1,91 1,70/1,79 1,79/1,85 2,07/2,29 2,22/2,47

Таблица 2. Теплофизические характеристики трубопровода

Обозначение, единица измерения Труба Изолятор Кожух

р, кг/м3 938,0 33,0 960,0

с, Дж/(кг -°С) 2300,0 1800,0 1880,0

Л, Вт/(м -°С) 0,35 0,04 0,42

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

Рис. 4. Приближение зависимости температур в трубопроводах от температуры воздуха: 1 — температура подающей трубы, 2 — температура обратной трубы, 3, 4 — аппроксимации.

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

Г1 (Га) = 67.592 - 1.1528 Га,

температуру обратного отопительного трубопровода определяем по эмпирической формуле

Г2(Га) = 42.964 - 0.3872 Га,

температура холодной воды постоянная Г3 = 12 0С. В отопительный период температура воды в трубопроводах определяется по вышеупомянутым эмпи-

рическим формулам, а в летний период задается фиксированная температура Т1 = Т2 = 21 0С. Переход с отопительного в летний период и обратно происходит при повышении и понижении средней температуры воздуха пяти последних дней до 8 0 С [19].

Температура воздуха была взята как среднесуточная из архива климатических данных. Суммарный коэффициент теплообмена воздуха с поверхностью грунта меняется со временем и зависит от месяца, начиная с января были использованы следующие значения: аа = 3.0, 3.0, 5.0, 8.0,14.0,14.0,14.0,14.0,10.0, 7.0, 5.0, 4.0.

5. Расчеты

Для исследования сходимости решения рассмотрены различные сетки по пространству и времени. Моделируется 36 суток, на 16-е сутки отопительный сезон заканчивается, т. е. температура воды в трубопроводах становится постоянной.

Рис. 5. Изменение погрешности по времени для разных временных шагов: 1 — т = 48 часов, 2 — т = 24 часов, 3 — т = 12 часов, 4 — т = 6 часовю

При исследовании зависимости погрешности от шага по времени для минимизации влияния погрешности от пространственной аппроксимации использовалась сетка эталонного решения. Шаг по времени принимал следующие значения: т = 48, 24, 12, 6 часов (рис. 5). При уменьшении шага по времени численное решение сходится.

Для исследования погрешности от пространственной аппроксимации шаг по времени взяли т = 3 часа. Рассматривались четыре варианта сеток: очень грубая (775 узлов, 1448 ячеек), грубая (2623 узлов, 5054 ячеек), средняя (8445 узлов, 16521 ячеек) и мелкая (32101 узлов, 63467 ячеек). На рис. 6 представлена динамика погрешности при различных вычислительных сетках. Погрешность уменьшается примерно на порядок при улучшении сетки в два раза. На

Рис. 6. Изменение погрешности по времени для различных сеток: 1 — 775 узлов, 2 — 2623 узлов, 3 — 8445 узлов, 4 — 32101 узловю

Рис. 7. Температура в скважине в различные даты по: 1 — натурным наблюдениям, 2 — результатам моделированияю

16-е сутки расчета отключается отопление и повышается погрешность решения около трубопроводов. Численное решение сходится при уменьшении стеки по пространству и времени.

Для численных экспериментов использовалась грубая сетка (2623 узлов, 5054 ячеек) с шагом по времени т = 12 часов. Численные расчеты сравнили с экспериментальными данными температур, полученными из скважины, расположенной на расстоянии 0.8 м от системы трубопроводов. Для сравнения использовались температуры четырех различных измерений в различные даты: 02.11.10, 05.05.12, 26.12.12, 17.05.13. Измерения температур первой даты

05.05.12

-10 -5 0 5 10

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

были взяты как начальные условия. Для всех остальных дат дается сравнение распределения температуры по глубине скважины (рис. 7). Распределение температуры по всей области для различных дат отображено на рис. 8.

Максимальное различие экспериментальных и расчетных данных получается для первого метра скважины. Для остального диапазона глубин скважины данные эксперимента и расчетов расходятся не более чем на 1 0 С.

ЛИТЕРАТУРА

1. Акимов М. П., Мордовской С. Д., Старостин Н. П. Численный алгоритм для исследования влияния бесканального подземного трубопровода теплоснабжения на вечномерзлые грунты // Мат. заметки ЯГУ. 2010. Т. 17, № 2. С. 125-131.

2. Акимов М. П. Математическое моделирование взаимодействия полимерного теплопровода теплоснабжения с вечномерзлыми грунтами // Научное творчество молодежи: Материалы XV Всерос. науч.-практ. конф. Ч. 1. Томск: Изд-во Том. ун-та, 2011. С. 258-260.

3. Васильев В. И., Максимов А. М., Петров Е. Е., Цыпкин Г. Г. Тепломассоперенос в промерзающих и протаивающих грунтах. М.: Наука, 1996.

4. Слепцов В. И., Мордовской С. Д., Изаксон В. Ю. Математическое моделирование тепло-обменных процессов в многолетнемерзлых горных породах. Новосибирск: Наука, 1996.

5. Васильев В. И., Попов В. В. Численное решение задачи промерзания грунта // Мат. моделирование. 2008. Т. 20, № 7. С. 119-128.

6. Васильева М. В., Павлова Н. В. Конечно-элементная реализация задачи замораживания фильтрующих грунтов // Мат. заметки ЯГУ. 2013. Т. 20, № 1. С. 202-212.

7. Павлов А. В. Расчет и регулирование мерзлотного режима почвы. Новосибирск: Наука, 1980.

8. Павлов А. В., Перльштейн Г. З., Типенко Г. С. Актуальные аспекты моделирования и прогноза термического состояния криолитозоны в условиях меняющегося климата // Криосфера Земли. 2010. Т. 14, № 1. С. 3-12.

9. Вабищевич П. Н. Численные методы решения задач со свободной границей. М.: Изд-во Моск. ун-та, 1987.

10. Самарский А. А. Введение в теорию разностных схем. М.: Наука, 1971.

11. Самарский А. А. Теория разностных схем. М.: Наука, 1977.

12. Стренг Г., Фикс Дж. Теория метода конечных элементов. М.: Мир, 1977.

13. Норри Д., де Фриз Ж. Введение в метод конечных элементов. М.: Мир, 1981.

14. Самарский А. А., Вабищевич П. Н. Вычислительная теплопередача. М.: Эдиториал УРСС, 2003.

15. Афанасьева Н. М., Васильева М. В., Захаров П. Е. Параллельное численное моделирование процесса заводнения нефтяного месторождения // Мат. заметки ЯГУ. 2011. Т. 18, № 1. С. 159-172.

16. Gaspar F., Grigoriev A., Vabishchevich P. Explicit-implicit splitting schemes for some systems of evolutionary equations // Int. J. Numerical Anal. Modeling. 2014. V. 11, N 2. P. 346-357. в 16 прошу указать год

17. Программа GMSH http://geuz.org/gmsh/.

18. Библиотека FEniCS http://fenicsproject.org/.

19. СНиП 41-02-2003. Тепловые сети / СНиП 41-02-2003. М.: Госстрой России, 2004.

Статья поступила 10 декабря 2014 г.

Акимов Мир Петрович, Захаров Петр Егорович Северо-Восточный федеральный университет им. М. К. Аммосова, Институт математики и информатики,

ул. Кулаковского, 48, Якутск 677000, Республика Саха (Якутия) mir_akimov@mail. ru

Захаров Петр Егорович

Северо-Восточный федеральный университет им. М. К. Аммосова, Институт математики и информатики,

ул. Кулаковского, 48, Якутск 677000, Республика Саха (Якутия) [email protected]

Матвеева Ольга Иннокентьевна

ОАО «Якутский государственный проектный

научно-исследовательский институт строительства»

ул. Дзержинского, 20, Якутск 677000, Республика Саха (Якутия)

[email protected]