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

ТЕХНИЧЕСКИЕ НАУКИ

УДК 519.6 В. С. Борисов

ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧИ О ПРОЦЕССАХ ПРОМЕРЗАНИЯ И ПРОТАИВАНИЯ В МНОГОЛЕТНЕМЕРЗЛЫХ ГРУНТАХ

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

Ключевые слова: уравнение теплопроводности, задача Стефана, вечная мерзлота, фазовые переходы, метод сквозного счета, методы конечных элементов, gnuplot, ParaView, FEniCS, GMSH.

V. S. Borisov

Numerical Solution of the Freezing and Thawing Processes Problem in Permafrost

Buildings in the permafrost are constructed taking into account the processes of thawing/freezing in the soil. Under the influence of heating season frozen ground under the building starts to melt, so foundation settlement could lead to the accident. Mathematical modeling of heat distribution in permafrost soils with phase transitions is described by the heat equation and Stefan condition. This paper presents the numerical implementation of this model in three-dimensional domain. In the formulation the model is written in the form of level set method without the explicit separation of the thawing front. Discretization is performed using the finite element method. The calculations were performed using the FEniCSlibrary for solver and GMSH program for mesh generation. A comparative analysis of the margin of error for the numerical and analytical solutions for one-dimensional formulation with Dirichlet and Neumann boundary conditions is made. The comparison is made on a different number of nodes, time step and smearing interval. Implemented computational algorithm can be applied with minimal changes for solving multidimensional problem. Geometric model for the problem in three-dimensional statement is constructed and results of numerical simulation for that problem are collected.

Keywords: heat equation, Stefan problem, permafrost, phase transitions, level set method, finite element methods, gnuplot, ParaView, FEniCS, GMSH.

БОРИСОВ Виктор Светославович - н. с. кафедры вычислительных технологий ИМИ СВФУ им. М.К. Аммосова.

E-mail: refactorman@gmail.com

BORISOV Victor Svetoslavovich - research Associate of Department of Computational Technologies, Institue of Mathematics and Information Science, North-Eastern Federal University named after M. K. Ammosov.

E-mail: refactorman@gmail. com

Введение

Здания на территории с вечной мерзлотой строятся с учетом процессов протаивания/промерзания в грунтах. Под влиянием теплоты сезонного отопления мерзлый грунт под зданием начинает подтаивать - происходит осадка фундамента здания, что может привести к аварийным последствиям. В условиях Крайнего Севера фундамент здания строится на базе свай -

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

Математическое моделирование процесса тепломассопереноса было описано в работах [1-5]. Типичной моделью с фазовыми переходами является задача Стефана. Для численного решения был разработан метод ловли фронта фазового перехода в узел разностной сетки. Постановка задачи Стефана с использованием единого уравнения теплопроводности для разных фаз с разрывными коэффициентами позволяет использовать метод сквозного счета [6]. В такой постановке задачу можно решить, используя как метод конечных разностей, так и метод конечных элементов [7-8].

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

Процесс распространения тепла с фазовыми переходами описывается следующим уравнением теплопроводности [3]:

ди

(сР + р1,8(и -ит))— = div(Xgrad и) + (x,г)

дг (1)

х € П, 0< г < Т,

где и = и (хД) - искомая функция распределения температурного поля, ит - температура фазового перехода, с - удельная теплоемкость, р = р (х) -плотность среды, L - энтальпия фазового перехода, д (и -и т) - дельта-функция Дирака, X = X (хД) - коэффициент теплопроводности, _Дх,0 - внутренние источники тепла, Т - конечное время. Геометрия О может быть представлена отрезком, прямоугольником или параллелепипедом в разных размерностях уравнения.

Коэффициенты в уравнении (1) - это кусочно-гладкие функции с разрывом в точке фазового перехода:

с5, и < ит ,

X —

К , и > ит , Х, и < ит ,

где с1=с1(х), Х1=Х1(х) - коэффициенты для талой зоны; с=сх(х), Хх=Хх(х) - для мерзлой зоны. Уравнение на границе Г = {Г1 и Г2 и Г3} дополняется смешанными граничными условиями первого, второго и третьего родов:

и (х,г ) = х,г), х € Г1,

ди

л — (х,г) = g2(х,г), х € Г2,

дп

ди

Л — (х, г) = а (и + gз(х, г)), х € Г3 дп

(2)

V =^ е Н (О): = а},

V ={\к е Н (О): w |г =0}.

(4)

Умножим уравнение (1), где и е V, на тестовую функцию V е у и проинтегрируем по области О с использованием формулы Грина:

I (ср + р1^8(и—ит))^Uvdx = ип дt

= — I (Xgrad и, grad v)dx +

и П

ди р ди р

X — vdx + I X — vdx + I IV ск.

Г дп о Г дп ¿П

ди

дп

дп

Интеграл по границе Г1 равен нулю по определению (4) пространства V. В интегралы по границам Г2, Г3 подставляем соответствующие значения из граничных условий (2). Таким образом, задача (1)-(3) ставится в вариационной постановке - найти функцию и е V, такую, что:

Г (ср + р^(и — ит))д- уйх = ¿я дt

= — Г (grad u,grad + ^ g2vdx + (5)

+Г а (и + g3)vdx + Г /V dx, V €

ди

Частную производной по времени — аппроксимируем неявной разностной схемой. Введем равномерную сетку по времени:

ю=(Г= тп, п= 0, 1,..., N1^ Т}.

Фазовый переход происходит в некой окрестности температуры фазового перехода. Для перехода к дискретному аналогу заменим дельта-функцию д ( и - и т) на дельтаобразную д (и - ит, Д)>0, где Д - величина полуинтервала размазывания [6]:

8(и - ит , Д) =

0, | и - ит |> Д,

2Д , 1 и - ит \< Д.

Разрывный коэффициент с заменим, используя интервал [ит- Д, ит+Д] для сглаживания значений. Здесь приведено сглаживание с применением линейной интерполяции:

где а - коэффициент теплообмена с внешней средой.

Распределение температуры в начальный момент времени задается следующим условием:

и(х,0) = и0(х), х е & и Г. (3)

Конечно-элементная дискретизация

Численное решение поставленной задачи (1)-(3) вычислим с помощью метода конечных элементов. Пусть Щ(О) - пространство Соболева, состоящее из функций V таких, что V2 и | Vw |2 имеют конечный интеграл в О. Определим V, V как подпространства Щ(О) следующим образом:

с, +

с, и < ит -А,

С - С, )(и - ит + А)

2А

^т \<А,

и > ит + А.

Коэффициенты р и X аппроксимируются аналогично.

Пусть Vh с V и у с у - сеточные конечномерные подпространства. Вариационная задача (5) приводится к следующей дискретной задаче. Найти сеточную функцию ип+1 € V из задачи:

с —

с

А

a(u"+\ v) = L(v), v eV h, n = 0,1,..., N-1 0

(6)

St,

r d^ du

pL — = —al — dt dx

- — dx

t >0

-1/2 .

du

XL— (0, t) = —q, t >0, q >0

dx

ft'

a,

= A_

PCL

as

= ЛS

Pcs

CL (uL - Um )

L :

StS

CS (um - US )

L :

и = и0 (х)

где формы а( , ) и £(•) заданы таким образом:

а(и,у) = Г с(и") — vdx + Г (Адgrad u,grad v)dx — Г auvdx,

^ П т П 3

— П

L(у) = Г с (и") — vdx + Г g2vdx + ^ g3vdx + ^ fvdx, ) = СдРд + Рд¿<5 (— — —т , Д).

Заметим, что функции f, g2, g3 приведены в записи только для удобства. В вычислениях используются их сеточные аналоги.

Задача в одномерной постановке

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

ди „ д2и „ „, ч

^ Р^Т = 0< х < ^)> г >0,

дt дх

ди . д2и ....

csР^Г = ЛУТТГ' £(0< х' t >0. дt дх

На подвижной границе £(0 установим температуру фазового перехода

и(Ш0= ит, г>0.

Поставим условие Стефана

ди

V =

Для задачи с граничным условием Дирихле точное решение имеет следующий вид

erf (

М( X ? - (ul - um )-

2Л laLt

)

erfc(

erf (k)

x

2Л laSt

0< x < %(t), t >0,

x > ), t >0,

u(X t) = us + {um - us )-

erfc(kv)

где число к находится из трансцендентного уравнения

Str

StS

— = к4л.

ek erf (к) ve k erfc(kv)

Для задачи с граничным условием Неймана точное решение [11] находят в таком виде:

К x 0 = —

erf (к)- erf x

2^1aLt

naL

0< x < l(t), t >0:

x

erfc

u(x,t) = us + (um -us)-

2yJ ast

x > l;(t), t >0,

и начальное значение температуры и положения границы

и(х,0)= и< ит, х> 0, С(0)= 0.

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

и(0,г)=и> ит, г>0. (7)

Второе - граничное условие Неймана, введем тепловой поток, пропорциональный г

erfc (kv)

где число к находится из трансцендентного уравнения:

q 1П 1 ^ =

рЧ aL ek ve 2 erfc(vk)

(8)

Приведем аналитическое решение одномерной задачи [9]-[10]. Движение фронта фазового перехода:

¡(г) = гк^а}, t >о.

Определим следующие константы для записи аналитических решений:

Для нахождения корней трансцендентных уравнений можно использовать итерационный метод Ньютона.

Численное решение одномерной задачи

Проведем численное решение задач в одномерной постановке, используя метод конечных элементов. Так, расчетная область ^ - это конечный отрезок 0< х < I. В то же время для задач в одномерной постановке область определения х - это полубесконечная прямая. Аппроксимируем бесконечную область, используя условие отсутствия потока:

ди

х—(1,0 = 0, t >0.

дх

В такой постановке численное и точное решение можно сравнивать до момента времени г<Ттхх Используя численное решение, вычислим это время экспериментально. Тогда при I = 10 м будет выполняться такое условие:

|ии(/,0-и(/,0|<10-6, 0<г<Ттхх= 23 суток.

Для численного решения возьмем следующие значения коэффициентов из [2]: р= 1400 кг/м3, с= 1710

x

)

Рис. 1. Зависимость погрешности от количества узлов: 1-16, 2-32, 3-64, 4-128, 5-256, 6-512

Рис. 2. Зависимость погрешности от шага по времени (с): 1-3600, 2-7200, 3-14400, 4-28800, 5-43200, 6-86400

Дж/(кг-°С), с = 1130 Дж/(кг-°С), Х1= 0,99 Вт/(м^°С),

1,33 Вт/(м-°С), 1=33,5 кДж/кг, и =2 °С, и=-5 °С, д=20411 Вт/м2, Д= 1 °С, /=10 м, Т=22 суток.

Для анализа численных решений определим следующие погрешности. Расчет абсолютной погрешности в пространстве 12:

АВЗ(;) = .

В анализе используем относительную погрешность в процентах:

АВЭО)

REL(t) =

^ и(х, ^)2 dx

:-100%.

Реализация вычислительного алгоритма сделана на основе библиотеки FEniCS [12], которая позволяет решать дифференциальные уравнения с частными производными методом конечных элементов. Из полученных решений строятся графики с использованием библиотеки gnuplot.

Проанализируем численное решение одномерной

задачи с граничным условием Дирихле. Рассмотрим поведение погрешности в зависимости от количества узлов сетки. На рис. 1 погрешность уменьшается при измельчении расчетной сетки. При количестве узлов более 128 погрешность не превышает 1,54 %. То есть измельчение сетки влечет за собой уменьшение погрешности.

Несмотря на чисто неявную схему, на погрешность сильно влияет выбор шага по времени. Расчеты с различными т проведены на сетке 512, Д=0,5 °С. На рис. 2 при слишком мелком т вычисляется больше временных шагов, что влечет за собой накопление ошибки. При большом т тоже увеличивается ошибка аппроксимации по времени. Изменение параметра Д и выбор шага по времени ограничили относительную погрешность до 0,93 %. При т= 14400 с накопленная ошибка на 0,1 % меньше, чем при т= 3600 с.

Параметр аппроксимации дельта-функции Д оказывает влияние на погрешность вычислений в зависимости от размера сетки. Для сетки с 512 узлами с шагом по времени т=14400 с. На рис. 3 представлена погрешность от параметра Д. Видно, что при разных значениях параметра сильно изменяется погрешность вычислений. Так,

ДНИ

Рис. 3. Зависимость погрешности от параметра Д (°С): 1-0,125, 2-0,25, 3-0,5, 4-0,75, 5-1,0, 6-1,25

Рис. 4. Сравнение погрешностей задач: 1 - Неймана, 2 - Дирихле

при Д=0,125 °С достигается минимальная погрешность в 0,16 %, но имеется скачок в графике. Следующий Д=0,25 °С выдает погрешность 0,45 %.

Решим одномерную задачу с граничными условиями Неймана, используя параметры из предыдущего анализа. Так, количество узлов =512, т= 14400 с, Д=0,25 °С Для таких параметров численное решение задачи Неймана выдает погрешность не более 2 %. Только для наглядности на рис. 4 представлено сравнение погрешности для задач Дирихле и Неймана с одинаковыми параметрами. Различие в граничном условии в точке х= 0 м. Погрешность вычислений хуже, чем для задачи Дирихле. Это предположительно обусловлено отсутствием известного значения температуры хотя бы в одной точке.

Численное решение трехмерной задачи

Проведем численное моделирование трехмерной задачи распространения тепла в грунтах с учетом фазовых переходов. Построим геометрию (рис. 5) с размерами 70x50x20 м3. На поверхности 2=0 м построим фундаменты соседних зданий с площадью 10x10 м2 и 14x9 м2. Расстояние между зданиями 5 м. Геометрия и расчетная сетка строится с помощью программы GMSH. Заметим локальное сгущение сетки (рис. 6) на поверхности г = 0 м и около фундаментов зданий.

Уравнение (1) дополняется граничными условиями:

и (х, t) = D, X е Г1, 0 < t < Т,

ди

X— (х, 0 = а(и — иаГ), х е Г2, 0<t < Т, дп

ди

X— (х,0 = 0, х е Г3, 0< t < Т, дп

ди

X — (х, t) = XG, х е Г4, 0< t < Т, дп

где постоянная температура, поддерживаемая в фундаменте здания, ^=15, коэффициент теплообмена а= 14, геотермический градиент 0= 0,027 °С/м и формула распределения температуры

ишг =-11- 35

+ 90)

365-24-3600

°С.

Возьмем следующее начальное условие: и0(х)= -5 °С.

Обработка результатов трехмерной задачи выполняется посредством программы для визуализации ParaView. На рис. 7 представлено распределение температуры через 10 лет на срезе х= 28 м. Белая линия - фронт фазового перехода.

Рис. 7. Температура срез по х= 28 м. через 10 лет. 1 - зима, 2 - весна, 3 - лето, 4 - осень

Заключение

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

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

Работа выполнена при финансовой поддержке гранта РФФИ №13-01-00719.

Л и т е р а т у р а

1. Carslaw H. S., Jaeger J. C. Heat in solids. - Clarendon Press, Oxford, 1959. - Vol. 1.

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

3. Samarskii A. A., Vabishchevich P. N. Computational Heat Transfer. - Wiley, 1996. - Vol. 1. - P. 418.

4. Cheng Q., Sun Y., Jones S. B., Vasilyev V. I., Popov V. V., Wang G., Zheng L. In situ measured and simulated seasonal freeze-thaw cycle: A 2-year comparative study between layered and homogeneous field soil profiles // Journal of Hydrology. - 2014. -Vol. 519. - P. 1466-1473.

5. Васильев В. И., Максимов А. М., Петров Е. Е., Цыпкин Г. Г. Математическая модель замерзания-таяния засоленного мерзлого грунта // Прикладная механика и техническая физика.

- 1995. - Т. 36, № 5. - С. 57-66.

6. Самарский А. А., Моисеенко Б. Д. Экономичная схема сквозного счета для многомерной задачи Стефана // Журнал вычислительной математики и математической физики. - 1965.

- Т. 5, № 5. - С. 816-827.

7. Вабищевич П. Н., Васильева М. В., Горнов В. Ф., Павлова Н. В. Математическое моделирование искусственного замораживания грунтов // Вычислительные технологии. - 2014.

- Т. 19, № 4. - С. 19-31.

8. Вабищевич П. Н., Васильева М. В., Павлова Н. В. Численное моделирование термостабилизации фильтрующих грунтов // Математическое моделирование. - 2014. - Т. 26, № 9.

- С. 111-125.

9. Tikhonov A. N., Samarskii A. A. Equations of Mathematical Physics. - Dover Publications, 2011. - P. 800.

10. Alexiades V, Solomon A. D. Mathematical modeling of melting and freezing processes. - Hemisphere Publishing Corporation, 1993.

11. Lozano R. F., Lara M. A. About the exact solution in two phase-stefan problem // Engenharia Termica (Thermal Engineering).

- 2007. - Vol. 6, no. 2. - P. 70-65.

12. Wells G., Mardal K. A., Logg, A. Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book. - Springer, 2012.

R e f e r e n c e s

1. Carslaw H. S., Jaeger J. C. Heat in solids. - Clarendon Press, Oxford, 1959. - Vol. 1.

2. Vasil'ev V. I., Maksimov A. M., Petrov E. E., Tsypkin G. G. Teplomassoperenos v promerzaiushchikh i protaivaiushchikh gruntakh. - M.: Nauka, Fizmatlit, 1996. - S. 224.

3. Samarskii A. A., Vabishchevich P. N. Computational Heat Transfer. - Wiley, 1996. - Vol. 1. - P. 418.

4. Cheng Q., Sun Y., Jones S. B., Vasilyev V. I., Popov V. V., Wang G., Zheng L. In situ measured and simulated seasonal freeze-thaw cycle: A 2-year comparative study between layered and homogeneous field soil profiles // Journal of Hydrology. - 2014. -Vol. 519. - P. 1466-1473.

5. Vasil'ev V. I., Maksimov A. M., Petrov E. E., Tsypkin G. G. Matematicheskaia model' zamerzaniia-taianiia zasolennogo merzlogo grunta // Prikladnaia mekhanika i tekhnicheskaia fizika.

- 1995. - T. 36, № 5. - S. 57-66.

6. Samarskii A. A., Moiseenko B. D. Ekonomichnaia skhema skvoznogo scheta dlia mnogomernoi zadachi Stefana // Zhurnal vychislitel'noi matematiki i matematicheskoi fiziki. - 1965. - T. 5, № 5. - S. 816-827.

7. Vabishchevich P. N., Vasil'eva M. V., Gornov V. F., Pavlova N. V. Matematicheskoe modelirovanie iskusstvennogo zamorazhivaniia gruntov // Vychislitel'nye tekhnologii. - 2014. - T. 19, № 4. - S. 19-31.

8. Vabishchevich P. N., Vasil'eva M. V., Pavlova N. V. Chislennoe modelirovanie termostabilizatsii fil'truiushchikh gruntov // Matematicheskoe modelirovanie. - 2014. - T. 26, № 9.

- S. 111-125.

9. Tikhonov A. N., Samarskii A. A. Equations of Mathematical Physics. - Dover Publications, 2011. - P. 800.

10. Alexiades V., Solomon A. D. Mathematical modeling of melting and freezing processes. - Hemisphere Publishing Corporation, 1993.

11. Lozano R. F., Lara M. A. About the exact solution in two phase-stefan problem // Engenharia Termica (Thermal Engineering).

- 2007. - Vol. 6, no. 2. - P. 70-65.

12. Wells G., Mardal K. A., Logg, A. Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book. - Springer, 2012.