Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt Моделирование температурных полей в вертикальной скважине с техногенной трещиной с использованием адаптивных сеток\ast Ю.О. Бобренева1, И.М. Губайдуллин1,2, Р.В. Жалнин3, В.Ф. Масягин3 Институт нефтехимии и катализа РАН1, Уфимский государственный нефтяной технический университет2, Мордовский государственный университет им. Н.П. Огарёва3 Предлагается методика на основе метода конечных объемов на локально адап- тируемых сетках для решения задачи о распределении температуры в нефтяном пласте с трещиной и вертикальной нагнетательной скважиной. Вычислитель- ный алгоритм реализован с использованием параллельной библиотеки BoxLib, предоставляющей структуры данных и подпрограммы для организации работы с адаптивными сетками. Данная библиотека поддерживает распараллеливание посредством MPI и OpenMP. В работе представлены результаты численного мо- делирования. Ключевые слова: вертикальная скважина, гидравлический разрыв пласта, тем- пературное поле, уравнение теплопроводности, адаптивные сетки, параллельные вычисления 1. Введение Развитие нефтяной и газовой промышленности России в последние десятилетия проис- ходит на фоне заметного ухудшения структуры запасов нефти и газа, что в основном свя- зано со значительной выработкой многих уникальных и крупных месторождений, а также вводом в разработку месторождений с трудноизвлекаемыми запасами. Степень выработ- ки запасов существенно зависит не только от правильного регулирования разработки с целью максимального извлечения остаточных запасов углеводородов, но и от полноты и достоверности информации о пласте и скважине. Одним из важнейших источников ин- формации являются гидродинамические (промысловые) исследования пластов и скважин. Совершенствование систем разработки нефтяных месторождений связано с применяемыми на промыслах мероприятиями по интенсификации добычи нефти. Промысловые исследо- вания скважин и пластов, поэтому приобретают все более важное значение как инструмент для оценки эффективности применяемых мероприятий. В процессе эксплуатации пластов и скважин исследования ведутся, главным образом, гидродинамическими методами, при этом уточняются характеристики пластов, выявляется эффективность мероприятий по воздей- ствию на призабойную зону пласта. Для более полного обеспечения информативности гид- родинамических методов исследований скважины, особенно для скважин с гидравлическим разрывом пласта [1], необходимо совместное рассмотрение гидродинамического состояния системы "скважина - пласт"с температурным полем, или термометрией. Термометрия ис- торически является первым методом исследования скважин. Основным параметром, кото- рый несет информационную нагрузку в данном методе, является температура. Температура – это энергетический параметр системы, и поэтому любое изменение системы вследствие изменения режима работы скважины, уменьшения или увеличения давления, промывки, нарушения целостности колонны и т.п. приводит к изменению температуры (распределе- \ast Работа выполнена при частичной финансовой поддержке РФФИ (Гранты № 15-07-01764 «Оптимальное управление химическими реакциями металлокомплексного катализа» и № 14-01-31260 мол_а «Исследование неявной методики для решения задач аэродинамики с использованием метода Галёркина с разрывными базисными функциями») 454 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt ния температуры) в скважине. Система скважина-пласт в этом отношении является очень чувствительной системой, т.к. на практике используются термометры с высокой разреша- ющей способностью. Профиль температуры зависит от скорости фильтрации, градиента давления, а также свойств жидкости и породы [2]. Для детального изучения термогидро- динамических процессов в пласте со скважиной необходима разработка математической модели совместной работы системы "скважина – трещина – пласт". Существующие высо- копроизводительные многопроцессорные компьютеры позволяют проводить многомерные расчеты по сложным физическим моделям, при этом, однако, требуемая пространственная дискретизация задачи не всегда может быть достигнута традиционными средствами [3]. Для более эффективного перераспределения вычислительных ресурсов существует подход, основанный на использовании блочных адаптивных сеток [4, 5], позволяющий контроли- ровать пространственный и временной шаг в заданной области решения. По сравнению с методами на неструктурированных сетках [6] или методами адаптивного измельчения, в ко- торых сетки всех уровней дискретизации принадлежат одному временному слою, и общий временной шаг задачи соответствует самой подробной сетке, данный подход позволяет ис- пользовать для каждого уровня детализации свой временной шаг, определяемый условием устойчивости численной схемы. Такая процедура адаптивного измельчения и перестраи- вания сеток, которая значительно экономит оперативную память и существенно ускоряет время моделирования по сравнению с равномерными сетками при том же пространственном шаге, реализована в библиотеке BoxLib [7, 8]. В настоящей статье впервые проведена моди- фикация и адаптация этой библиотеки для моделирования температурных полей в пласте с техногенной трещиной и вертикальной нагнетательной скважиной. 2. Метод локальной адаптации сеток В основе алгоритма адаптивного измельчения сеток (AMR - Adaptive Mesh Refinement) лежит метод Бергера и Олигера [4] построения иерархии вложенных сеток, первоначально предложенный в 1984 г. для задач, описываемых системой дифференциальных уравнений гиперболического типа. Под адаптивным построением вложенных сеток понимают рекурсивный процесс из- мельчения ячеек самой грубой сетки до тех пор, пока на уровне с самой мелкой сеткой решение не будет найдено с заданной точностью, причем измельчение происходит только в тех областях, в которых это необходимо (большие градиенты, разрывы, контактные грани- цы). При моделировании нестационарных задач на сетках разных уровней дискретизации периодически оценивается ошибка вычислений и, если критерий ошибки в некоторых ячей- ках перестает удовлетворяться, то эти ячейки помечаются как “требующие измельчения” и добавляются в сетки более высокого уровня измельчения. Наоборот, некоторые ячейки мелкой сетки удаляются из нее, если высокое разрешение в данной области больше не тре- буется. Таким образом, иерархия сеток автоматически адаптируется в ходе интегрирова- ния уравнений так, чтобы области с большими градиентами параметров получали большее разрешение по времени и пространству. Принципиальным недостатком метода AMR явля- ется его относительная сложность по сравнению с методами, использующими одну сетку. Так, для реализации AMR требуются рекурсивные алгоритмы и иерархические структу- ры данных. Поэтому для решения физических задач целесообразно использовать готовые библиотеки [7], которые берут на себя основную вычислительно-алгоритмическую работу по поддержанию иерархии вложенных сеток, межпроцессорному обмену при параллельных вычислениях и т.д. В данной работе использовалась библиотека BoxLib [8], предоставляю- щая функциональность для использования AMR и разработанная группой ученых из Center for Computational Sciences and Engineering (CCSE) at Lawrence Berkeley National Laboratory, California, USA. Явный вид граничных условий не принципиален для нижеизложенного алгоритма. Ре- 455 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt шение задачи находится рекурсивно путем прохода по сеткам всех уровней l от самого грубого l = 0 до самого мелкого l = lmax . Если ячейка уровня l покрыта ячейками более мелкого уровня l + 1, то значения величин в такой ячейке заменяются усредненными зна- чениями с уровня l + 1. На границе сеток, принадлежащих разным уровням, производится коррекция потоков для грубых ячеек, примыкающих к более мелкой сетке. Шаги по вре- мени для уровней l и l + 1 связаны соотношением \Delta tl+1 = nl1 \Delta tl , где nlref – параметр ref измельчения сетки уровня l по отношению к l + 1. 3. Конечно-объемная дискретизация уравнения теплопровод- ности Рассматривается математическая модель, описывающая распределение температурного поля в системе «скважина-трещина-пласт» в процессе эксплуатации вертикальной нагнета- тельной скважины [9]. Температура закачиваемой жидкости меньше температуры пласта. Распределение температурного поля в техногенной трещине и продуктивном пласте описывается уравнениями: \biggl( \biggr) \partial Tf \partial \partial Tf \partial Tf \alpha f t = \lambda f t - cl \rho l Vf , t > 0, 0 < x < Xf , (1) \partial t \partial x \partial x \partial x \biggl( \biggr) \biggl( \biggr) \biggl( \biggr) \partial Tm \partial \partial Tm \partial \partial Tm \partial Tm \partial Tm \alpha mt = \lambda mt + \lambda mt - cl \rho l Vmx + Vmy , (2) \partial t \partial x \partial x \partial y \partial y \partial x \partial y t > 0, 0 < x < Lx , 0 < y < Ly , с начальными и граничными условиями: Tm | t=0 = Tf | t=0 = T0 , Tf | x=0 = T0 + \Delta T, Tm | x=Lx = Tm | y=Ly = T0 , \partial Tm \partial Tm \partial x | x=0 = 0, \partial y | y=0 = 0, \partial Tf \partial Tf \lambda mt \partial T \partial Tm \partial x | x=Xf = \lambda f t \partial x | x=Xf , \lambda mt \partial y | y=Wf /2 = \lambda f t \partial y | y=Wf /2 , m Tm | x=Xf = Tf | x=Xf , Tm | y=Wf /2 = Tf | y=Wf /2 , где \alpha f t – объемная теплоемкость в трещине, \alpha mt – объемная теплоемкость в пласте, \lambda f t – коэффициент теплопроводности в трещине, \lambda mt – коэффициент теплопроводности в пласте, cl – удельная теплоемкость жидкости, \rho l – плотность жидкости, Vf – скорость фильтрации жидкости в трещине, Vmx , Vmy – скорость фильтрации жидкости в пласте, Tf – температура в трещине, Tm – температура в пласте, T0 – температура в начальный момент времени. Для дискретизации уравнений (1)–(2) на равномерной прямоугольной сетке исполь- зуется противопоточная схема первого порядка точности для конвективных членов [10] и центральные разности для диффузионных членов: d(Tf )i (Tf )i - 1 - 2(Tf )i + (Tf )i+1 (T\^f )i+1/2 - (T\^f )i - 1/2 \alpha f t = \lambda f t - cl \rho l V f , i = 1, \=Nf , (3) dt \Delta x2 \Delta x d(Tf )ij (Tf )i - 1j - 2(Tf )ij +(Tf )i+1j (T ) - 2(T ) +(T ) \alpha mt = \lambda mt \Delta x2 + \lambda mt f ij - 1 \Delta yf 2ij f ij+1 - (4) dt \biggl( \biggr) (T\^f )i+1/2j - (T\^f )i - 1/2j (T\^ ) - (T\^ ) - cl \rho l Vmx \Delta x + Vmy f ij+1/2\Delta y f ij - 1/2 , i = 1, \=Nx , j = 1, \=Ny , где (T\^f )i+1/2j , (T\^f )ij+1/2 – численные конвективные потоки через границу между ячейками с индексами (i, j) и (i + 1, j) и ячейками с индексами (i, j) и (i, j + 1), соответственно. 456 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt Также вводится в рассмотрение слой фиктивных ячеек, значения в которых задаются исходя из граничных условий. Для дискретизации по времени полученной системы дифференциально разностных уравнений (3) – (4) используется явная схема первого порядка точности. 4. Параллелизм библиотеки BoxLib Распараллеливание алгоритмов выполнено за счет встроенных средств параллельной библиотеки BoxLib . Данная библиотека поддерживает гибридную модель на базе MPI и OpenMP. BoxLib автоматически выполняет декомпозицию расчетной сетки на каждом уровне l на доступное количество процессоров. Также библиотека предоставляет набор под- программ для формирования набора фиктивных ячеек и выполнения межпроцессорного обмена для их заполнения значениями с соседних процессоров. Разработчики гарантируют возможность масштабирования более чем на 200000 процессоров. Основной абстракцией, реализующей параллелизм, является тип данных MultiFab [8]. Это структура, инкапсулирующая набор массивов данных FAB, соответствующих каждо- му блоку сетки. Любая применяемая к данным типа MultiFab операция выполняется для каждого массива FAB согласно их распределению между процессорами. Данные MultiFab на каждом уровне адаптации распределены по процессорам независимо. Библиотека под- держивает параллельный ввод/вывод и параллельный рестарт. 5. Результаты численного моделирования Рис. 1. Схематичное изображение системы «скважина-трещина-пласт» В начальный момент времени, когда скважина еще не запущена в эксплуатацию, во всей области задаются параметры: начальная температура T0 =363.15 К, удельная теплоемкость c1 =2000 Дж/(кг\cdot К), плотность жидкости \rho 1 950 кг/м3, объемная теплоемкость в трещине \alpha f t =1.316676\cdot 106 Дж/(м3\cdot K), объемная теплоемкость в пласте \alpha f t =1.316676\cdot 106 Дж/(м3\cdot K), коэффициент теплопроводности в трещине \lambda f t =2.5208, коэффициент теплопроводности в пласте \lambda mt =2.5208. Длина пласта 2\cdot Lx =300 м, ширина пласта 2\cdot Ly =50 м, толщина трещины Wf =0.005 м, полудлина трещины Xf =40 м. На правой и верхней границах задаются следующие температуры T1 = Th = 363.15 К, скорость фильтрации жидкости сонаправлена с радиус-вектором точки. В каждой ячейке скорость постоянна и равна значению в ее геометрическом центре. В начале координат за- 457 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt дается температура Tf | x=0,y=0 = 293.15 К. На рис. 1 представлено схематичное изображение исследуемой системы. На рис. 2 представлена сетка, построенная с помощью библиотеки BoxLib, на начальный момент времени. Рис. 2. Фрагмент расчетной сетки на начальный момент времени t = 0 Рис. 3. Распределение температуры на момент времени t = 4.0 На рис. 3 представлено распространение температурного фронта по пласту. Отметим, что ввиду симметрии рассматриваемой области и граничных условий достаточно рассмат- ривать только четверть расчетной области. Из рисунка видно, как во время работы на- гнетательной скважины холодная закачиваемая жидкость охлаждает пласт. Значительные изменения температуры наблюдаются вблизи скважины и вдоль распространения трещины, в частности на створках трещины. На рис. 4 – 5 представлены фрагменты расчетной сетки и уровни адаптации на раз- личные моменты времени. Подробная сетка формируется только в области с высоким гра- диентом температуры, что позволяет экономить оперативную память и сократить время расчета. 458 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt Рис. 4. Фрагмент расчетной сетки на момент времени t = 0.054 Рис. 5. Фрагмент расчетной сетки и уровни адаптации на момент времени t = 2.321 459 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt 6. Заключение Проведенные вычислительные эксперименты показали, что разработанная компьютер- ная программа позволяет проводить моделирование температурных полей в системе «сква- жина - трещина - пласт». Реализованная в библиотеке BoxLib процедура адаптивного из- мельчения сеток показала свою эффективность и стабильность, в расчетах использовано до 5 уровней измельчения. Наряду с использованием параллельных вычислительных техноло- гий, использование адаптивных сеток дает заметный выигрыш по времени выполнения и использованию оперативной памяти по сравнению с равномерными сетками с шагом, соот- ветствующим шагу самого подробного уровня. Литература 1. Малышев А.Г. и другие. Анализ влияния технологических факторов и механических свойств горных пород на эффективность ГРП. // В. кн. «Нефть Сургута». -М.: Нефтяное хозяйство, 1997. С. 224-237. 2. Чекалюк Э.Б. Термодинамика нефтяного пласта.-Москва: «Недра»-1965, 238 с. 3. Губайдуллин И.М., Линд Ю.Б., Коледина К.Ф. Методология распараллеливания при решении многопараметрических обратных задач химической кинетики // Вычислительные методы и программирование: новые вычислительные технологии. 2012. Т. 13. № 2 (26). С. 28 – 36. 4. Berger M.J., Oliger J. Adaptive mesh refinement for hyperbolic partial differential equations // J. of Computational Physics. 1984. 53. 484–512 5. Поварицын М.Е., Захаренков А.С., Левашов П.Р., Хищщенко К.В. Моделирование многокомпонентных гидродинамических течений с использованием адаптивных сеток // Вычислительные методы и программирование: новые вычислительные технологии. 2012. Т. 13. № 3. С. 424 – 433 6. Масягин В.Ф., Жалнин Р.В., Тишкин В.Ф. О применении разрывного конечно-элементного метода галëркина для решения двумерных уравнений диффузионного типа на неструктурированных сетках // Журнал Среневолжского математического общества. 2013. Т. 15. № 2. С. 59 – 65. 7. Dubey A., Almgren A., Bell J., Berzins M., Brandt S., Bryan G., Colella P., Graves D., Lijewski M., Loffler F., O’Shea B., Schnetter E., Van Straalen B., Weide K.. A Survey of High Level Frameworks in Block-Structured Adaptive Mesh Refinement Packages // Journal of Parallel and Distributed Computing, 74, pp. 3217-3227, 2014 8. Bell J., Almgren A., Beckner V., et al. BoxLib User’s Guide. URL: https://ccse.lbl.gov/BoxLib/BoxLibUsersGuide.pdf (дата обращения: 15.02.2016). 9. Рамазанов А.Ш., Нагимов В.М. Аналитическая модель для расчета температурного поля в нефтяном пласте при нестационарном притоке жидкости // Электронный научный журнал «Нефтегазовое дело». 2007. №1. С.1-9. 10. Куликовский А.Г., Погорелов Н.В., Семенов А.Ю., Математические вопросы численного решения гиперболических систем уравнений, изд. 2-е, исправл. и доп., Физматлит, М., 2012, 656 с. 460 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt Modelling of temperature fields in a vertical well with a man-made crack using adaptive grids Yu.O. Bobreneva1, I.M. Gubaydullin1,2, R.V. Zhalnin3, V.F. Masyagin3 Institute of petrochemistry and catalysis of the Russian Academy of Science1, Ufa State Petroleum Technological University2, Ogarev Mordovia State University3 FVM-based method using AMR is proposed for solving the problem about temperature distribution in the oil reservoir with crack and vertical wells. The computational algorithm uses the BoxLib parallel computation library, that contains data structures and subroutines for adaptive meshes managing. The results of numerical modelling are presented in the work. Keywords: vertical wells, hydraulic fracturing, temperature field, heat conduction equation, AMR, parallel computing References 1. Malyshev A.G. and etc. Analiz vliyaniya tekhnologicheskikh faktorov i mekhanicheskikh svoystv gornykh porod na effektivnost’ GRP [Analysis of the impact of technological factors and mechanical properties of rocks on the efficiency of hydraulic fracturing]. // In book ”Neft surguta” [”Surgut oil”]. Moscow. Publishing of the ”Neftyanoe xozyajstvo”, 1997. p. 224 – 237. 2. Chekalyuk E.B. Termodinamika neftyanogo plasta [Thermodynamics oil reservoir]. Moscow, Publishing of the ”Nedra”, 1965. 238 p. 3. Gubajdullin I.M., Lind Yu.B., Koledina K.F. Metodologiya rasparallelivaniya pri reshenii mnogoparametricheskikh obratnykh zadach khimicheskoy kinetiki [Parallelization methodology for solving multiparameter inverse problems of chemical kinetics] // Vychislitelnye metody i programmirovanie: novye vychislitelnye texnologii [Numerical methods and programming: new computing technologies]. 2012. V. 13. No. 2 (26). P. 28 – 36. 4. Berger M.J., Oliger J. Adaptive mesh refinement for hyperbolic partial differential equations // J. of Computational Physics. 1984. 53. 484–512 5. Povaritsyn M.E., Zakharenkov A.S., Levashov P.R., Khishchshchenko K.V. Modelirovanie mnogokomponentnykh gidrodinamicheskikh techeniy s ispol’zovaniem adaptivnykh setok [Simulation of multi-hydrodynamic flows using adaptive mesh refinement] // Vychislitel’nye metody i programmirovanie: novye vychislitel’nye tekhnologii [Numerical methods and programming: new computing technologies]. 2012. V. 13. No. 3. P. 424 – 433 6. Masyagin V.F., Zhalnin R.V., Tishkin V.F. O primenenii razryvnogo konechno-elementnogo metoda galerkina dlya resheniya dvumernykh uravneniy diffuzionnogo tipa na nestrukturirovannykh setkakh [On the application of the discontinuous Galerkin finite element method for solving two-dimensional equations of diffusion type on unstructured meshes] // Zhurnal Srenevolzhskogo matematicheskogo obshchestva [Journal of the Middle-Volga mathematical society]. 2013. V. 15. No. 2. p. 59 – 65. 7. Dubey A., Almgren A., Bell J., Berzins M., Brandt S., Bryan G., Colella P., Graves D., Lijewski M., Loffler F., O’Shea B., Schnetter E., Van Straalen B., Weide K.. A Survey of High Level Frameworks in Block-Structured Adaptive Mesh Refinement Packages // Journal of Parallel and Distributed Computing, 74, pp. 3217-3227, 2014 461 Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016) agora.guru.ru/pavt 8. Bell J., Almgren A., Beckner V., et al. BoxLib User’s Guide. URL: https://ccse.lbl.gov/BoxLib/BoxLibUsersGuide.pdf (accessed: 15.02.2016). 9. Ramazanov A.Sh., Nagimov V.M. Analiticheskaya model’ dlya rascheta temperaturnogo polya v neftyanom plaste pri nestatsionarnom pritoke zhidkosti [Analytical model to calculate the temperature field in the oil reservoir under unsteady fluid inflow] // Elektronnyj nauchnyj zhurnal «Neftegazovoe delo» [Electronic scientific journal ”Oil and Gas Business”]. 2007. No. 1. p.1-9. 10. Kulikovskij A.G., Pogorelov N.V., Semenov A.Yu. Matematicheskie voprosy chislennogo resheniya giperbolicheskikh sistem uravneniy [Mathematical problems in the numerical solution of hyperbolic systems], ed. 2nd, corrected. and ext., Moscow. Publishing of ”FIZMATLIT”, 2012. 656 p. 462