ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ 3D ЗАДАЧ ТЕПЛОПРОВОДНОСТИ С ФАЗОВЫМИ ПЕРЕХОДАМИ НА ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМАХ С РАСПРЕДЕЛЕННОЙ ПАМЯТЬЮ Текст научной статьи по специальности «Физика»

УДК 536.2:519.87:622.276

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ 3D ЗАДАЧ ТЕПЛОПРОВОДНОСТИ С ФАЗОВЫМИ ПЕРЕХОДАМИ НА ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМАХ С РАСПРЕДЕЛЕННОЙ ПАМЯТЬЮ

И. В. Бычин1, Т. В. Гавриленко2, В. А. Галкин3, А. В. Гореликов4, А. В. Ряховский5

Сургутский государственный университет, 1 igor-bychin@yandex.ru, 2 taras.gavrilenko@gmail.com, 3 val-gal@yandex.ru, 4 gorelikov_a@list.ru, 5 echo47@rambler.ru

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

Ключевые слова: сопряженный теплообмен, фазовый переход, MPI.

NUMERICAL SIMULATION OF 3D HEAT CONDUCTION PROBLEMS WITH PHASE TRANSITIONS USING AN HPC CLUSTER

I. V. Bychin1, T. V. Gavrilenko2, V. A. Galkin3, A. V. Gorelikov4, A. V. Ryakhovsky5

Surgut State University, 1 igor-bychin@yandex.ru, 2 taras.gavrilenko@gmail.com, 3 val-gal@yandex.ru,

4 gorelikov_a@list.ru, 5 echo47@rambler.ru

The article describes the iterative algorithm for solution of conjugate heat conduction problems with phase transitions. The phase transition boundary is determined by using explicit method (accurate to a control volume), the temperature field in every phase is computed by implicit method. The control volume method in the Cartesian coordinate system and the enthalpy method were used. The algorithm is implemented in software for numerical solution of freezing-thawing problem in layered heterogeneous media. The software is parallelized using MPI standard. Tests were carried out and demonstrated the agreement between results, analytical solutions and numerical solutions obtained by other researchers.

Keywords: conjugate heat transfer, phase transition, MPI.

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

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

Математическую модель процессов промерзания-растепления слоисто-неоднородного массива пород составляет:

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

2) двухфазная задача Стефана.

Далее представлена математическая постановка двухфазной задачи Стефана. В «жидкой» фазе (талый грунт):

дТ

РЪ-^ = Т) + Б. (1)

В «твёрдой» фазе (мерзлый грунт):

Т

РСз-Т^ = ёгаб Т) . (2)

На границе раздела фаз температура равна температуре фазового перехода:

Т = Т (3)

и справедливо уравнение баланса энергии для поверхности раздела фаз (условие Стефана):

(4)

где А — удельная теплота плавления; #п — нормальная к поверхности раздела фаз составляющая скорости движения данной поверхности, п — нормаль к поверхности раздела фаз, внешняя по отношению к твердой фазе.

Далее представлена энтальпийная формулировка задачи. Объемная энтальпия вводится следующим образом:

Н(Т) = [ С(Т) йТ + Ьв(Т-Т), (5)

■)Т0

\ о,е< о '

где С объемная теплоемкость, Ь — объемная теплота фазового перехода, Т0 — произвольно выбранное значение температуры. Уравнение теплопроводности в каждой фазе записывается в виде:

^ = div(kgraйT)+ Б. (6)

о1

По известной энтальпии, температура может быть найдена по формуле:

Н <Н2 : Т = НСН2 + Т Т = \ Н2 < Н < Нх : Т = Т . (7)

Н>НХ: Т = + Т

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

1. Дискретный аналог уравнения энтальпии (6) решается по явной схеме:

Н(п) = Ф (г(п~ ^ Ы + Н(п~, (8)

где А I = 1п — 1 - шаг по времени; Н(п) - сеточное значение энтальпии на текущем временном слое ¿п; Ф (Т(п~- конечно-разностный оператор аппроксимирующий диффузионные потоки на гранях контрольного объема и объемные источники тепла в уравнении (6). Коэффициент теплопроводности на гранях контрольных объемов рассчитывается, как среднегармоническое значение между двумя ближайшими значениями в узлах сетки.

2. Первое приближение на временном слое 1п для поля температуры Т* определяется из соотношения:

Н <Н2: Т = НсН1 + т Г = 4 Н2 < Н < Н1 : Т = Т . (9)

Н > Н1 : Т = НСН1 + Тг

Граница фазового перехода определяется с точностью до контрольного объема. Контрольный объем, температура которого равна Т* (или энтальпия которого лежит в интервале H < H < H ) считается двухфазным. Теплоемкость и коэффициент теплопроводности в двухфазном объеме определяется как средневзвешенное значение между соответствующими значениями в твердой и жидкой фазах в зависимости от доли жидкой фазы.

3. Следующее приближение для температуры Т** находится в результате итерационного решения дискретного аналога уравнения теплопроводности для каждой из фаз (1), (2). В данном случае используется полностью неявная схема для дискретизации по времени:

т** _ т(п~ 1)

С-—-= Ф(Г*) . (10)

Уравнение (10) совместно с уравнениями состояния С(Т), k(T) решаются до достижения численной сходимости на данном временном слое. После достижения сходимости, Т(n) = Т**.

Таким образом в разработанном алгоритме по явной схеме определяется только граница фазового перехода (с точностью до контрольного объема), а поле температуры в каждой фазе рассчитывается по полностью неявной схеме. Такой подход позволяет ослабить ограничения на шаг по времени по сравнению с явной схемой и организовать пересчет на данном временном слое всех теплофизических характеристик С(Т), k(J).

Программный код комплекса распараллелен с помощью функций стандарта MPI [4, 5]. При распараллеливании использовался метод двухмерной декомпозиции по координатам x, y. Расчетная область разбивалась на перекрывающиеся подобласти в каждой из которых задача решается отдельным процессом. Каждый процесс пересылает значения из импортируемых контрольных объемов и получает значения в теневых контрольных объемах от соседних процессов (рис. 1).

ООО

Oil

neo

Рис. 1. Двухмерная декомпозиция расчетной области (вид сверху на сечение расчетной области плоскостью xOy). Черным цветом обозначены теневые грани, синим цветом импортируемые грани, серым - грани на которых заданы граничные условия.

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

Расчеты с использованием MPI-версии программного комплекса проводились на аппаратно-программном комплексе АПК-5 СурГУ (пиковая производительность 5 Тфлопс, 16 вычислительных узлов, конфигурация узла: 2 x Intel XEON CPU E5-2650 v2, 32 Гб ОЗУ. Общее количество физических ядер: 256). Предварительные расчеты на сетке 512x512x182 контрольных объемов с двухмерной де-

Рис. 2. Блок-схема вычислительного ядра программного комплекса.

композицией расчетной области на 256 подобластей демонстрирует ускорение в 80 раз по сравнению с однопоточным кодом.

Программный комплекс протестирован на задачах с автомодельными решениями [6]. Проведено сравнение с результатами расчетов [6], полученных методом с явным выделением подвижных границ. На задачах с автомодельными решениями итерационная схема, реализованная в программном комплексе, демонстрирует второй порядок аппроксимации по пространственным переменным.

Тест 1. Одномерная однофазная задача Стефана с аналитическим решением

Рассматривается задача о плавлении длинного тонкого полубесконечного бруска 0 < х < оо, расположенного вдоль оси Ох, с адиабатической боковой поверхностью. На левом конце задается тепловой поток. Задача однофазная, в том смысле, что уравнение теплопроводности решается только для жидкой фазы, температура твердой фазы в течении всего процесса плавления равна температуре фазового перехода. $(£)-положение границы раздела фаз вдоль оси Ох. В жидкой фазе 0 < х < $(1) :

Т _ д2Т

д1 _ дх2'

Условия Стефана на границе раздела фаз х _ $(0 :

(11)

А--)

Ч 9х):

В «твёрдой» фазе $(0 < х

Начальные условия £ _ 0:

Граничные условия х _ 0:

Аналитическое решение:

Т _ 0, ( ) _ 1.

Т = 0.

Т _ 0 ; $(0) _ 0.

Т

дх

_ — е.

Зт(0 _ Тап(х, 0 _

{

0 < х < $ап(£) : е*-х - 1 $апУ) < х : 0

(12)

(13)

(14)

(15)

(16)

Далее представлены результаты численного решения данной задачи, на сетках 20 и 40 контрольных объемов вдоль оси Ох.

На данном тесте разработанная схема демонстрирует второй порядок аппроксимации по пространственным переменным (рис. 3-4).

Тест 2. Протаивание слоисто-неоднородного массива пород вокруг скважины Далее представлена задача о протаивании слоисто-неоднородного массива пород вокруг скважины [6]. В пределах каждого пласта перенос тепла описывается двухфазной задачей Стефана. В зоне протаивания:

С1 (Т )&Т _ Шу (к1 (Т) Т) + ^. (17)

В мерзлой породе:

На границе протаивания:

д1

(Т)|Т _й\у(к$(Т) рай Т) .

Т _ Т(,

\дп)

тН _ ^п

(18)

(19)

I

где 1?п нормальная к поверхности раздела фаз составляющая скорости движения данной поверхности, п нормаль к поверхности раздела фаз, внешняя по отношению к твердой фазе Ц — теплота фазового перехода на единицу объема породы данного слоя, определяется соотношением: Ц _ ^(Щ — Щ)) Ь, где Щ- общая весовая влажность, Щ — содержание незамерзшей воды, 7 — вес скелета породы, Ь — удельная теплота фазового перехода лёд-вода Ь _ 334 • 103Дж/кг.

а) б)

Рис. 3. Распределение температуры едоль стержня е момент еремени t = 0,25 (красный цеет — аналитическое решение, зелёный — численное): а) результаты для сетки 20 контрольных объемое; б) результаты для сетки 40 контрольных объемое

а) б)

Рис. 4. Распределение температуры едоль стержня е момент еремени t = 0,85 (красный цеет — аналитическое решение, зелёный — численное): а) результаты для сетки 20 контрольных объемое; б) результаты для сетки 40 контрольных объемое

На контактах пластов задаются условия сопряжения:

[T ] =0,

■kd_r

dz

= 0.

(20)

(21)

На нижней и верхней границах массива задаются адиабатические граничные условия:

T

z = 0, z = h, — = 0.

z

В начальный момент времени t = 0 массив заполнен мерзлой породой при температуре Tin (Tin < Tf). В массиве находятся изотермические тепловые источники (скважины), температура в стволе скважины T (T > Tf). Температура газа в стволе скважины задавалась равной T0 = 40°C, температура протаива-ния пород Tf = 0°C. Внешний диаметр скважины — 499 мм, диаметр эксплуатационной колоны — 168

а) б)

Рис. 5. Распределение температуры на безразмерное время t = 0,85: а) аналитическое решение; б) численное решение

мм. Теплопроводность заполняющего межколонные пространства цемента принимались равной 0,35 Вт/(м • °C). Теплофизические свойства пород приведены в табл. 1. Строение разрезов, принятых в качестве исходных данных, соответствует реальным характеристикам мерзлой толщи на Бованенковском газоконденсатном месторождении. Начальная температура пород принималась равной Tin = —4,5°C. Моделирование проводилось на периоды 5-летней и 30-летней эксплуатации газодобывающих скважин. Расчетная сетка: nx = 256, ny = 256, nz = 182.

Таблица 1

Исходные характеристики пород, приняты в расчетах

Тип пород , г/см3 W, % Wo, % Л, Вт/(м • °C) C, кДж/(м3 • °C)

талые мерзлые талые мерзлые

лед - - - 0.57 2.22 4200 2090

песок 1.45 28 0.9 1.5 1.8 2570 1830

глина 1.35 32 0.8 1.35 1.65 2260 1970

На рис. 6 представлены результаты решения задачи растепления двухслойного массива пород вокруг скважины за 30 лет. Верхняя часть массива мощностью 15 метров состоит из чистого льда, а нижняя из мерзлой глины мощностью 30 метров (т.е. высота расчетной области 45 м). Размеры расчетной области вдоль осей Ох и Оу составляют 26.6 м. Задача обладает осевой симметрией относительно оси скважины, это позволяет (при расчетах в прямоугольных координатах) расположить скважину вдоль вертикального ребра расчетной области. На вертикальных гранях расчетной области, прилегающих к скважине, из-за осевой симметрии задачи задаются адиабатические граничные условия.

На рис. 7 представлены результаты решения задачи растепления трехслойного массива пород вокруг газодобывающей скважины за 5 лет. Верхняя часть массива мощностью 11м состоит из глины, нижняя из песка мощностью 30 м, между ними пласт льда мощностью 4 м.

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

Далее на рис. 8-9 представлены примеры, демонстрирующие возможности программного комплекса по решению задач растепления от нескольких тепловых источников и по решению задач растепления от наклонных скважин.

Рис. 6. Решение задачи растепления двухслойного массиеа пород еокруг скеажины за 30 лет: а) изопоеерхность температуры Т=0 (фронт плаеления) и изотермы е сечениях по оси 02; б) изотермы е сечении по оси Оу

а) б)

Рис. 7. Решение задачи растепления трехслойного массиеа пород еокруг скеажины за 5 лет: а) изопоеерхность температуры Т=0 (фронт плаеления) и изотермы е сечениях по оси 02; б) изотермы е сечении по оси 0у

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

Работа выполнена при поддержке РФФИ, гранты: № 13-01-12051 офи_м, № 15-41-00013 р_урал_а.

а) б)

Рис. 8. Пример решения задачи растепления трехслойного массива пород от наклонной скважины за 1 год: а) изоповерхность температуры Т=0 (фронт плавления) и изотермы в сечениях по оси О2; б) изотермы в сечении по оси Оу

Рис. 9. Пример решения задачи растепления трехслойного массива пород от трех скважин (изоповерхность температуры Т=0 и изотермы в сечениях по оси О2): а) полгода; б) 1 год; в) 5 лет; г) 30 лет

ЛИТЕРАТУРА

1. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М. : Энерго-атомиздат, 1984. 152 с.

2. Versteeg H. An Introduction to Computational Fluid Dynamics: The Finite Volume Method. Harlow : Prentice Hall, 2007. 520 p.

3. Pletcher R. H. Computational Fluid Mechanics and Heat Transfer. Boca Raton : CRC Press, 2012. 774 p.

4. Gropp W. Using MPI: portable parallel programming with the message-passing interface. Cambridge : MIT Press, 2014. 336 p.

5. Антонов A. С. Технологии параллельного программирования MPI и OpenMP. М. : Издательство Московского университета, 2012. 344 с.

6. Сигунов Ю. А. Методы решения классической задачи Стефана. Сургут : РИО СурГПУ, 2009. 138 с.