О численном моделировании теплопередачи от скважины к пласту породы в пакете Ansys

Размещено на http://www.allbest.ru/

О ЧИСЛЕННОМ МОДЕЛИРОВАНИИ ТЕПЛОПЕРЕДАЧИ ОТ СКВАЖИНЫ К ПЛАСТУ ПОРОДЫ В ПАКЕТЕ ANSYS

Автор(ы): Меньшакова Татьяна Вячеславовна магистрант, факультет математики и информационных технологий, Федеральное государственное бюджетное образовательное учреждение высшего образования «Национальный исследовательский Мордовский государственный университет им. Н. П. Огарёва»

Сыромясов Алексей Олегович кандидат физико-математических наук, доцент, кафедра прикладной математики, дифференциальных уравнений и теоретической механики, факультет математики и информационных технологий, Федеральное государственное бюджетное образовательное учреждение высшего образования «Национальный исследовательский Мордовский государственный университет им. Н. П. Огарёва»

Аннотация

В статье рассматривается простейшая задача о передаче тепла от скважины к окружающему ее пласту породы. Полученное ранее аналитическое решение задачи применяется для построения области, в которой далее проводятся вычисления с использованием пакета ANSYS Workbench. Оценивается суммарная погрешность и скорость сходимости численного метода.

1.Постановка и аналитическое решение задачи

Имеется пласт нефтесодержащей породы, в центре которого пробурена цилиндрическая скважина (рис. 1).

Рис. 1. Схематичное изображение пласта нефтесодержащей породы, имеющего скважину

В нее под давлением закачивается вода, температура которой постоянна. Поскольку она не равна начальной температуре пласта, между скважиной и пластом возникает тепловой поток. Требуется найти распределение температуры в пласте [1].

Рассмотрим задачу в простейшей безразмерной постановке: будем считать, что тепло распространяется лишь вдоль оси Ox^*. Тогда пласт задается неравенством x^*>0, а граница скважины ? равенством x^*=0. Описанная геометрия показана на рис. 2. Обозначим безразмерную температуру через T^*, время - через t^*.

Рис. 2. Схематичное изображение упрощенной постановки задачи

Распределение температуры в пласте породы описывается уравнением теплопроводности

которое имеет начальное условие:

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

Здесь безразмерный параметр Bi есть число Био, вычисляемое по формуле:

где в -- коэффициент теплоотдачи,

R -- радиус скважины,

л -- коэффициент теплопроводности в пласте.

К этим уравнениям необходимо присоединить условие затухания по времени на бесконечности:

В результате мы имеем систему:

Для решения задачи были взяты следующие конкретные данные: л = 2,5208 Вт/(м ? К), в = 20 … 100 Вт/(м²? К), R = 0,05 м. Далее нам также понадобятся плотность жидкости в пласте с = 950 кг/м³ и объемная теплоемкость пласта б = 1316676 Дж/ (м³ ? К).

Решение системы (1), полученное интегральным преобразованием Лапласа по времени, имеет вид [5]:

где используется дополнительная функция ошибок [2]:

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

2.Численное решение задачи в пакете ANSYS

Для численного решения задачи (1) используем систему ANSYSWorkbench, в которой выберем тип расчета TransientThermal.

Нами была создана геометрия задачи: мы нарисовали прямоугольник размерами x^* = 75 м и y^* = 1 м описали вещество с ранее указанными плотностью, теплопроводностью и удельной теплоёмкостью c_б = б/с. Вещество с этими характеристиками было приписано созданной геометрии.

На внешней границе было заданограничное условие типа Temperature: T^*= 1, внутри пласта было задано начальное условие T^*= 1. На левой границе пласта, т.е. на границе со скважиной, мы задали условие типа Convection. В качестве его параметра было взято число Био, принимающее значения от 0,396699 до 1,9835 при выбранных нами в, R, л. В самой скважине температура всегда равна нулю.

Для более точного решения было произведено сгущение сетки у границы скважины и пласта c помощью инструмента Inflation (рис. 3).

Рис. 3. Схематичное изображение геометрии задачи в пакете ANSYS

Посовещавшись с коллегами из Уфы, мы выяснили, что размерное время t, в течение которого нас интересует расчет, составляет 10 дней. Исходя из этого, мы взяли безразмерное время t^*= лt/(бR²), изменяющееся от 0 до 650.

Как только мы упростили задачу, заменив бесконечную область на конечную, наше решение приобрело неустранимую погрешность, связанную с неточностью исходных данных, принятых для численного расчета. Как показывает подстановка выбранных нами значений x^*и t^* в (2), эта погрешность должна составлять не более 3-4%. Такая погрешность вполне допустима для практических целей.

Как видно из рис. 4, с течением времени распределение температуры меняется, начиная с левого края области.

Рис. 4. Распределение температуры в конце расчета

Для получения таблицы значений температуры с шагом по времени t^* = 50 с по координате x^*= 7,5 м мы выводили результаты вдоль контура с помощью инструмента Path пакета ANSYS.

3.Сравнение аналитического и численного решений

Будем сравнивать результаты нашего численного решения с полученным ранее аналитическим решением (2). Нами рассмотрены три случая численной реализации задачи: шаг вычислений по координате составлял Дx^*=1; Дx^*=0,5; Дx^*=0,25. Используем число Фурье, характеризующее соотношение между скоростью изменения тепловых условий в окружающей среде и скоростью перестройки поля температуры внутри рассматриваемой системы; оно зависит от размеров тела (или шага вычислений) и коэффициента его температуропроводности. В нашем случае, когда уравнения обезразмерены, число Фурье Fo вычисляется по формуле

и называется также параболическим числом Куранта [3].

Изменяя Дx^*, мы должны изменять и Дt^* так, чтобы величина Fo оставалась постоянной. Отсюда Дt^*= Fo ? (Дx^*)², и при Fo = 5 для наших Дx^* шаг по времени Дt^* составляет соответственно 5; 1,25; 0,3125. Такое значение числа Фурье соответствует документации пакета ANSYS.

В каждом из случаев абсолютная погрешность в точке вычисляется как модуль разности температур, полученных при численном и аналитическом решениях. Ниже представлен фрагмент таблицы (табл. 1) при Дx^*=1, показывающий, что абсолютная погрешность не превышает 0,03?0,04; это соответствует относительной погрешности в 3-4%, о которой было сказано ранее.

скважина температура погрешность адаптивный

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

Далее оценим, насколько сильно наши решения в целом отличаются друг от друга. Для этого необходимо выбрать метрику (расстояние) в функциональном пространстве, которому принадлежат решения. Нами было проверено два вида метрик.

Пусть имеются функции f и g, заданные в области D. Тогда метрика может вычисляться одним из следующих способов [4]:

В нашем случае область D задается парами (t^* , x^*), вместо f и g будут взяты аналитическое и численное решения задачи (1).

Использовать (3) нам не удалось. Аналитическое решение (2) вычисляется точно и при x^* > ? меняется очень слабо, а при численном решении задачи условия на правой границе области заданы заранее и не изменяются с течением времени. Вместе данные обстоятельства приводят к тому, что максимум разности температур, полученных при численном и аналитическом решениях, не изменяется. Поэтому полученные с помощью (3) значения метрики не характеризуют точность решения в целом.

Применить (4) напрямую тоже не получается, потому что численное решение задачи известно только в некоторых точках области D. Поэтому заменим в (4) интеграл на сумму и в качестве суммарной погрешности ДT_total будем рассматривать выражение

где Д и Д_i,j - модуль разности между аналитическим и численным решениями в произвольной точке (t^* , x^*) и в точке (t^*_i , x^*_j), соответственно.

Указанная суммарная погрешность складывается из погрешности математической модели, неустранимой погрешности, погрешности вычислительного метода и погрешности округления. Известно, что погрешность вычислительного метода должна быть в 2-5 раз меньше неустранимой погрешности [6].

Нами найдено, что при Дx^*=1 величина (5) равна 0,1322549; при Дx^*=0,5 она составляет 0,0335669, а при Дx^*=0,25 она равна 0,0088048.

Эти данные позволяют оценить порядок сходимости численного решения. Говорят, что порядок сходимости равен n, если при шаге вычислений h

Тогда при шаге вычислений h/2

откуда

В качестве h будем рассматривать Дx^*. Сравнивая суммарные погрешности при h = 1 и h = 0,5, получим n ? 1,98; аналогично, при h = 0,5 и h = 0,25 получим n ? 1,93. Среднее значение n ? 1,95 слабо отличается от обеих найденных величин, т.е. примененная схема вычислений устойчива.

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

Список использованных источников

1. Бобренева Ю. О., Губайдуллин И. М., Жалнин Р. В., Масягин В. Ф. Моделирование температурных полей в вертикальной скважине с техногенной трещиной использованием адаптивных сеток // Параллельные вычислительные технологии (ПаВТ'2016): труды междунар. науч. конф. (Архангельск, 28 марта - 01 апреля 2016 г.). Архангельск, 2016. С. 454-462.

2. Лыков А.В. Теория теплопроводности. М.: Высш. шк., 1967. 600 с.

3. Аристова Е. Н., Лобанов А. И. Практические занятия по вычислительной математике в МФТИ: учеб. пособие. Ч. II. М.: МФТИ, 2015. 310 с.

4. Колмогоров А. Н., Фомин С. В. Элементы теории функций и функционального анализа. М.: Наука, 1976. 544 с.

5. Сыромясов А. О., Меньшакова Т. В. Об аналитических методах моделирования теплопередачи в системе “скважина + пласт” // Дифференциальные уравнения и их приложения в математическом моделировании: материалы XIII Междунар. науч. конф. (Саранск, 12-16 июля 2017 г.). Саранск: СВМО, 2017. С. 322-329.

6. Кафедра физики твердого тела Петрозаводского государственного университета [Электронный ресурс]: сайт. - URL: http://dssp.petrsu.ru/p/tutorial/meth_calc/files/02.shtml (дата обращения: 18.03.2018.).

Размещено на Allbest.ru