Конечно-элементная реализация задачи замораживания фильтрующих грунтов Текст научной статьи по специальности «Физика»

УДК 519.63:517.958

КОНЕЧНО-ЭЛЕМЕНТНАЯ РЕАЛИЗАЦИЯ ЗАДАЧИ ЗАМОРАЖИВАНИЯ ФИЛЬТРУЮЩИХ ГРУНТОВ

М, В, Васильева, Н, В, Павлова

1. Введение

Исследование изменений температурного режима грунтов является необходимым элементом инженерно-геологического обоснования строительства объектов в районах вечномерзлых грунтов. При сезонном оттаивании мерзлых грунтов изменяются их физико-механические свойства, что приводит к нарушению несущей способности фундаментов зданий и сооружений. Для укрепления фундаментов оснований зданий и сооружений используется метод замораживания грунтов. Замораживание грунтов производится с помощью специальных холодильных установок или с помощью сезонных охлаждающих устройств, не требующих затрат электрической энергии. Использование сезонных охлаждающих устройств также позволяет производить охлаждение грунтов в районах, где электричество недоступно, например, на нефтепроводах и газопроводах [1].

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

©2013 Васильева М. В., Павлова Н. В.

Рис. 1. Расчетная область.

[2-4] и в предыдущей работе авторов [5].

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

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

Рассмотрим математическую модель, описывающую распределение температуры при наличии фазовых переходов твердая фаза — жидкая фаза при некоторой заданной температуре фазового перехода Т* в области П = П- и (рис. 1). Здесь П+(£) — область, занятая жидкой фазой, где температура превышает температуру фазового перехода:

«+(*) = {х | х € П, Т(х,£) > Т*}, и П-(£) — область, занятая твердой фазой:

«-(*) = {х|х€П, Т(х, £) < Т * }.

Фазовый переход происходит на границе раздела фаз 5 = £(£).

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

поглощением и выделением скрытой теплоты [6]:

(а(ф) + р+Ьф')(^- - - сНу(А(</>) ёгас1 Т) = 0, (1)

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

Коэффициенты уравнения удовлетворяют следующим соотношениям:

а(ф) = р-с- + ф(р+с+ - р-с-), \(ф) = \~ + ф(Х+ - Х-),

О при Т<Т*,

Ф= ,

1 1 при Т > Т*,

где — плотность и удельная теплоемкость талой и мерз-

лой зоны соответственно.

Поскольку рассматривается процесс распространения тепла в пористой среде, имеем

с-р- = (1 - m)cscPscJr тСгРг, с+р+ = (1 - ш)с8ср8о + тсшрш,

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

X- = (1 - т)Х^ + тХI, Х+ = (1 - т)Х^ + тХш.

На практике фазовые превращения не происходят мгновенно и могут происходить в малом интервале температуры [Т* - А, Т* +Д]. В качестве функции ф можно взять фд:

0 при Т < Т* - Д, ФА = { Т~2А+Л при I Д • Г ■ /' — Д.

1 при Т > Т* + Д,

О при Т < Т* - Д,

Ф'А={ ш при Т*-А<Т<Т* + А, О при Т > Т* + Д.

Тогда получим уравнение для температуры в области Л

(дТ \

(а(фЛ) + Р1Ьф'А)( _+ugradTj - й1у(Х(фл) ёгас1Т) = 0, (2)

которое является стандартным нелинейным параболическим уравнением.

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

Т(х, о) = То, X € п, т = тс, X € тв, -к^- = о, х е г/гв. (з)

дп

Здесь Гд — место контакта с замораживающей колонкой.

Для учета фильтрации в грунте запишем уравнение для определения пластового давления в случае полностью насыщенной несжимаемой пористой среды:

— = х € П+, (4) где А = Уравнение (4) дополняется граничными условиями

к др к др —— = здг, х е Гдг, р(х) = зс, х 6 гв, —т- = о, х е 5. (5)

^дп ^дп

Здесь дП+ = 5и ГN иГ5 — подвижная граница фазового перехода.

Полученная задача (4), (5) является задачей с подвижной границей Для численного решения такой задачи без перестроения сетки воспользуемся методом фиктивных областей, который основывается на переходе к решению задачи в более широкой области [6,7]. Приближенное решение, зависящее от параметра продолжения е, будем искать во всей расчетной области П. При использовании варианта метода фиктивных областей с продолжением по старшим коэффициентам решение определяется из уравнения

— с1пг(Аб ^адр^ 0, х€П. (6) Здесь разрывный коэффициент Аб(х) определяется так:

АГ , Г А, X € П +, АЛх) = <

1 Ае, х € О

при достаточно малом е. Задание подобным образом коэффициентов имеет физический смысл и моделирует фильтрацию в области Л- с очень малым коэффициентом Ае ^ 0 при е ^0.

Для уравнения (6) формируются граничные условия, соответствующие граничным условиям исходной задачи (4), (5): к др

— — (х) = дм, X е Гдг, р(х) = до, X е Гд, Мдп (п

к дР ^ ' —/(х) = 0, хе 5П\(ГдгПГв).

М дп

3. Вариационная формулировка

Проведем дискретизацию полученной системы уравнений (2), (3) и (6), (7) с использованием метода конечных элементов [8]. Умножим уравнения для температуры и давления на тестовые функции уг> ур и проинтегрируем с использованием формулы Грина:

J(а(фА) + рЬф'А) +

п

+ У (А(фд) ёга<1 Т, ёга<1 Ут) ¿х = 0 УУТ € Н (П), (8) п

J(Ае ёгас1р, ёгас1 Ур) ¿х = 0 € Н(П). (9)

п

Здесь Н(П) — пространство Соболева, состоящее из функций У таких, что У2 и |УУ|2 имеют конечный интегр ал вАи Н(П) = {У € Н(П) | У|Гл = 0}.

Для аппроксимации по времени уравнения для температуры применим стандартную чисто неявную схему и воспользуемся простейшей линеаризацией (с предыдущего временного слоя):

/уп+1 _ тп /к

I (а(ф1)+р1Ьф'1){ТП+1т Г" + ^ёгас1р",ёгас1Г"+1^г;т^ п

+ J (А(Фд) ёга<1Тп+1 ,ёга<Ьт) ¿х = 0, (10)

J(\egradpn+1 ,gradvp) dx = 0. (11)

n

Для численного решения необходимо перейти от непрерывной вариационной задачи (10), (11) к дискретной задаче. Введем конечномерные пространства Vh, Vh, (Vh, Vh G H1(íl)) и определим в них дискретную вариационную задачу: найти Th,ph G Vh такие, что

J(а(ф1) + Р1Ьф'1) (^Тр1 + gradpi, grad dx

Q

+ J (А (Фд ) grad T"+1 , grad vT) dx n

= TJИФд) + Р'^Ф'Д)ThnvT dx VvT G Vh, (12) n

J (Ae gradp^1 , grad v£) dx = 0 Vv£ G Vh. (13)

n

Заметим, что выбор пространства Vh непосредственно вытекает из вида применяемых конечных элементов.

4. Вычислительная реализация

Процесс численного решения поставленной задачи состоит из следующих этапов:

• построение геометрии и генерация сетки;

•

•

В качестве модельной задачи рассмотрен процесс искусственного замораживания фильтрующих грунтов с использованием замораживающей колонки, на которой поддерживается постоянная температура T

T с

Геометрическая область строится с помощью программы Net gen, которая генерирует тетраэдальную сетку в трехмерной расчетной области. На рис. 2 изображена трехмерная расчетная область протяженностью 4 м по каждому направлению, в середине области находится

Рис. 2. Расчетная сетка.

замораживающая скважина радиусом 0.1 м. Расчетная сетка построена со сгущением вблизи замораживающей колонки и содержит 441 653 ячеек.

Для численного решения задачи используется программный пакет FEniCS, реализующий метод конечных элементов [8]. Полученные результаты, значения температуры и давления на каждом временном слое записываются в файл и визуализируются с помощью программы Par avie w.

5. Численный эксперимент

Приведем результаты вычислительного эксперимента. Расчеты проводились при входных параметрах задачи, представленных в табл. 1.

Результаты численных расчетов через 20 дней работы замораживающей колонки в трехмерной области представлены на рис. 3, 4. Расчет проводился при г = 1 день и Tmax = 20 дней.

На рис. 5 и б представлены изотермы при T = 0 и распределения температурного поля для различных моментов времени (5, 10, 20 дней) в срезе по направлениям y и z.

Численное исследование эффективности распараллеливания на вычислительном кластере «Ариан Кузмин» СВФУ приводится в табл. 2. Время счета на вычислительном кластере на 16 потоках составило около 30 секунд, при запуске на 4 потока — около 75 секунд, а на одном потоке — 260 секунд, что показывает достаточно хорошую эффективность распараллеливания вычислительного процесса. Декомпозиция области распараллеливания на 16 потоках представлена на рис. 7.

Таблица 1. Параметры задачи

Обозначение Значение Метрика Описание

Тс -30.0 град Температура замораживающей колонки

То 2.0 град Начальная температура

Т * 0.0 град Температура фазового перехода

Ро 1.0е6 Па Начальное давление

Р1 1.1е6 Па Давление на левой границе

Рг 1.0е6 Па Давление на правой границе

т 0.2 Пористость

к 1.0е-13 м2 Проницаемость грунта

М 1.0е-3 Па/сек Вязкость воды

Ь 1.04 * 1.0е8 Дж/кг Удельная теплота фазового перехода

Срэс 1.20*1.0е6 Дж/м3 Объемная теплоемкость грунта

cрf 4.20*1.0е6 Дж/м3 Объемная теплоемкость воды

Срг 1.91*1.0е6 Дж/м3 Объемная теплоемкость льда

Авс 1.1 Вт/(м град) Коэффициент теплопроводности грунта

А 0.56 Вт/(м град) Коэффициент теплопроводности воды

Аг 2.26 Вт/(м град) Коэффициент теплопроводности льда

Рис. 3. Распределение температуры через 20 дней.

Рис. 5. Сравнение температурного поля для различных моментов времени: 5, 10 и 20 дней (срез по у).

Рис. 6. Сравнение температурного поля для различных моментов времени: 5, 10 и 20 дней (срез по г).

Рис. 7. Декомпозиция области распараллеливания на 16 потоках.

Таблица 2. Зависимость времени счета от количества запущенных процессов

Количество процессов 1 2 4 8 16

Время счета (сек) 260.79 151.4 75.17 41.85 29.70

ЛИТЕРАТУРА

1. Мат. междунар. науч.-практ. конф. по инженерному мерзлотоведению, посвященной ХХ-летию создания ООО НПО «Фундаментстройаркос». Тюмень: Сити-Пресс, 2011.

2. Сапунов П. Е. Исследование процесса промерзания полусферической подземной ледопородной емкости, размещенной в фильтрующем пласте // Численные методы решения задач фильтрации многофазной несжимаемой жидкости. Новосибирск, 1980. С. 228-235.

3. Бондарев Э. А., Васильев В. И. Искусственное замораживание фильтрующих грунтов // Численные методы решения задач фильтрации многофазной несжимаемой жидкости. Новосибирск, 1987. С. 38-47.

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

5. Васильева М. В., Павлова П. В. Численное решение нестационарной задачи искусственного замораживания фильтрующих грунтов // Мат. заметки ЯГУ. 2010. Т. 17, №2. С. 113-123.

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

7. Вабищевич П. П. Метод фиктивных областей в задачах математической физики. М.: Изд-во Моск. ун-та, 1991.

8. Logg A., Mardal К.-А., Wells G. N. Automated solution of differential equations by the finite element method. The FEniCS book. Berlin: Springer-Verl., 2012. (Lect. Notes Comput. Sei. Eng.; V. 84).

г. Якутск

5 февраля 2013 г.