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

М.М. Хасанов, А.Н. Ситников, А.А. Пустовских, А.П. Рощектаев, Н.С. Исмагилов, Г.В. Падерин, Е.В. Шель, ООО «Газпромнефть НТЦ», Санкт-Петербург, Россия

Источник: Журнал «Георесурсы»

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

Как всякие большие системы, объекты нефтегазодобычи требуют использования целой иерархии моделей – от дифференциальных до интегральных, от детерминированных до адаптивных, способных описать не только различные уровни организации систем, но и взаимодействие между этими уровнями (Дьячук, 2003). решением всех этих задач по иерархичному моделированию процессов нефтегазовой отрасли занимается научный инжиниринг (англ. SciencEngineering) – область науки, находящаяся на стыке системного инжиниринга, фундаментальных наук и кибернетики (рис. 1).

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

Геологическое моделирование: применение спектрального подхода для моделирования геофизических полей

Разработка спектрального подхода к задачам геологического моделирования была впервые представлена в работе (Байков и др., 2010), в которой авторами предложен способ, основанный на разложении каротажа на коэффициенты по ортонормированному базису в пространстве функций, интегрируемых с квадратом L2. Дальнейшее развитие метод получил в работах (Байков и др., 2012; Хасанов и др., 2015), в которых осуществлена более детальная проработка теоретических основ спектрального метода, а также представлены некоторые результаты расчётов, основанных на спектральном методе геологического моделирования. стоит также упомянуть о том, что спектральный анализ каротажа получил применение для целей литолого-фациального анализа (Хасанов и др., 2014).


Рис. 1. Научный инжиниринг (SE) – на стыке областей различных походов

Спектральный метод геологического моделирования – метод моделирования трёхмерных кубов геофизических свойств на основе скважинных данных. Принятая в спектральном методе геологического моделирования математическая модель представлена моделируемой областью стохастическим полем G(x,y,h), определенным на этой области и заданным на некотором вероятностном пространстве. Для такой модели скважинные данные для вертикальных скважин с фиксированной координатой по латерали (x*, y*) представляют собой случайный процесс, параметризованный переменной, характеризующей глубину G(x*,y*,h) = f(h). если в моделировании участвуют  N скважин с координатами (x ,y ), на которых задано моделируемое свойство функциями fi(h), то эти функции можно рассматривать как известные значения некоторой реализации стохастического поля G(x,y,h).

Общий принцип спектрального метода заключается в последовательной реализации нескольких шагов: разложение функций fi(h) по некоторому ортонормированному базису, моделирование коэффициентов разложения в межскважинном пространстве, восстановление моделируемого стохастического поля по этим коэффициентам в каждой точке области D.

Известно, что в гильбертовом пространстве интегрируемых с квадратом функций (L2) существуют ортонормированные базисы, которые позволяют разложить любую функцию, заданную в этом пространстве, в ряд по коэффициентам разложения, определяемым однозначно. Функции, определяющие моделируемое свойство, как функции с конечной энергией принадлежат этому пространству. Выберем в пространстве L2 базис многочленов Лежандра Pj(h) который используем для разложения функций fi(h):

где коэффициенты c i определяются как скалярное произведение в пространстве L2. Набор коэффициентов разложения c i для каждого уровня разложения j определён в точках (xi,yi), которые являются координатами скважин, и является известными значениями некоторой реализации стохастического поля cj(xi,yi). 

Для моделирования коэффициентов c i для всех (x,y) используется спектральный метод моделирования (Пригарин, Михайлов, 2005), основанный на теореме об интегральном представлении случайных процессов и полей.

Построение обусловленной реализации смоделированного поля c (x ,y ) известными значениями c i осуществляет

В результате осуществления вышеприведённых шагов получается набор реализаций стохастических полей коэффициентов разложения cj(xi,yi), которые обусловлены скважинными значениями c i. Для восстановления реализации всего поля G(x,y,h) в каждой точке моделируемой области D осуществляется суммирование базисных функций по смоделированным коэффициентам:


В силу построения реализация поля является обусловленной скважинными данными и воспроизводящей статистические характеристики исходных данных.

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

Спектральный метод в силу алгоритма своего построения является принципиально параллелизуемым, что позволяют легко масштабировать вычисления. Кроме того, метод реализует построение модели без привязки к сетке (т.н. grid-free simulation), что позволяет осуществлять моделирование на сетках любой конфигурации и сложности, производить измельчение части сетки и уточнение модели на этой области и, наконец, параллельно генерировать реализации моделируемого свойства в различных областях модели. Прирост производительности спектрального моделирования на типовой рабочей станции геолога-модельера составляет от 2 до 4 раз по сравнению с классическими методами в зависимости от размера генерируемой модели. 

Применение спектрального метода на практике продемонстрировало его способность хорошо воспроизводить геофизические поля в межскважинном пространстве и в неразбуренных областях, в особенности при сравнении с традиционными методами моделирования (Рис. 2) (Хасанов и др., 2014; Пригарин, Михайлов, 2005).


Рис. 2. Сравнение реальной и синтетической кривых гамма каротажа. DGK, красная кривая, – реальный каротаж, SP_ DGK, синяя кривая, – синтетически каротаж, смоделированный спектральным методом, SGS_DGK, чёрная кривая, – синтетический каротаж, смоделированный методом SGS.

Гидродинамическое моделирование: Метод источников в оптимизационных технологических расчетах. 

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

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

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

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

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

В представленной работе рассмотрен пример использования метода источников/стоков для решения задачи выбора оптимального времени отработки нагнетательных скважин (Ситников и др., 2015). Время отработки – это период работы скважины в режиме добычи перед переводом ее в нагнетание. Распределение давления от вертикальной скважины, для которой задана динамика жидкости, записывается в виде интеграла Дюамеля:


где –  коэффициент пьезопроводности, P0 – начальное пластовое давление, µ – вязкость нефти, K – проницаемость, h – толщина пласта, С – общая сжимаемость, m – пористость, r – расстояние от источника, t – время расчета, τ – переменная интегрирования. Задача об определении дебита жидкости скважины, действующего при постоянном забойном давлении, сводится к отысканию функции из интегрального уравнения (3). Для кусочно-постоянного дебита жидкости выражение (3) трансформируется в уравнение (4):


где qn – кусочно-постоянный дебит на каждом шаге, Pc n – забойное давление, rc – радиус скважины, n – количество временных шагов, – интегральная показательная функция, Dt – шаг по времени. Отметим, что при выводе уравнения 4 использовалось решение линейного источника, которое, в свою очередь, получено в предположении малости радиуса ствола скважины. В подавляющем большинстве практических ситуаций такое приближение обосновано. 

Используя принцип суперпозиции (что обеспечивается линейностью уравнения пьезопроводности), можно рассчитать динамику нестационарных дебитов для системы из N скважин:


Для трещин ГРП конечной проводимости система уравнений (5) должна быть дополнена системой уравнений по расчету забойных давлений на источниках, имитирующих трещины. 

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


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

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

Достоверность расчетов предложенным методом проверялась сравнительным анализом динамики добычи жидкости с результатами расчетов в коммерческом программном продукте и аналитическими зависимостями, полученными в статье (Хасанов и др., 2013). Апробация выполнялась на основе данных о строении одного из месторождений компании ПАО «Газпромнефть». Полученная зависимость безразмерной накопленной дисконтированной добычи от времени отработки нагнетательной скважины показана на рис. 3. Видно, что наибольшая дисконтированная добыча может быть получена при условии использования нагнетательных скважин в добыче на протяжении примерно девяти месяцев. Более ранний или поздний перевод скважин окажется экономически менее выгодным. Данный результат не только подтвердился сериями расчетов на полномасштабной гидродинамической модели, но и находится в согласии с характерным временем перевода скважин, которое было получено опытным путем при эксплуатации месторождения (Ситников и др., 2015). 

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

Рис. 3. Зависимость безразмерной дисконтированной добычи Qdsc/max(Qdsc ) от времени перевода скважины в нагнетание T (Ситников и др., 2015) 

Геомеханика: сравнение параметров ГРП в безразмерных переменных по дизайну и гидродинамическим исследованиям 

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

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

Из графика, приведенного на рис. 4 видно, что явной зависимости между длиной трещины и объемом закачки в общем случае не наблюдается. Причины, которые к этому привели, достаточно очевидны, и заключаются в том, что помимо объёма закачки сшитого геля (и связанной с ним массы проппанта) на длину, очевидно, могут влиять такие параметры пласта, как модули Юнга по разрезу, мощности пластов, сжимающие напряжения, перпендикулярные трещине, коэффициенты трещиностойкости, а также технологические параметры дизайна ГРП – реология жидкости, скорость закачки и концентрация проппанта. 

Для анализа полученных данных требуется найти выборку скважин, которая была бы полностью идентична по всем параметрам ГРП кроме какого-то одного, после чего данную процедуру следовало бы повторить для всех параметров ГРП, чтобы установить степень влияния каждого из них. Определив данные степени влияния, можно было бы установить наиболее значимые параметры для геометрии трещины ГРП. Это бы позволило предположить параметры, ошибка в которых могла привести к значительному систематическому расхождению между длиной трещины, запланированной по дизайну, и длиной трещины, полученной по гидродинамическим исследованиям скважины. На практике, однако, подобрать такие выборки трещин ГРП, которые не отличались бы по всем параметрам кроме какого-то одного, не представляется возможным. 

Все кейсы выбраны так, чтобы закачиваемые жидкости ГРП были идентичны, в том числе по концентрации полимера. Тогда можно говорить о схожей реологии жидкостей ГРП для данных случаев. Таким образом, хотя бы по этому параметру данные по операции ГРП идентичны.


Рис. 4. Зависимость длины трещин по дизайну и ГДИС от объема закаченной жидкости

Для метода, предложенного в этой статье, это является необходимым условием. Скорости закачки жидкости ГРП при дизайне также изменяются слабо, так что на большой разброс и образование «облака» на рис. 4 этот фактор не мог оказать решающее воздействие. 

Таким образом, остаются в основном разнообразные геомеханические параметры, которые могут меняться вдоль пласта достаточно значительно (особенно мощность пласта). Для успешного статистического анализа требуется метод, который понизит размерность задачи по этим параметрам, и приведёт трещины ГРП в пластах с разными геомеханическими свойствами к одному «знаменателю». Таким методом является введение безразмерных параметров задачи о развитии трещины ГРП. 

Развитие трещины по модели Planar3D в случае закачки неньютоновской жидкости описывается тремя законами (Хасанов и др., 2017): 

- законом Гука; 

- законом вязкого трения; 

- законом сохранения массы. 

Обезразмеривание уравнений проведем аналогично тому, как описано в работе (José I. Adachi et al., 2010) для случая трещины без проппанта; также будем считать, что утечки в пласт отсутствуют, модуль Юнга – однородный, а литология – трехслойная и симметричная. В результате получаем следующий список безразмерных параметров:Рис. 4. Зависимость длины трещин по дизайну и ГДИС от объема закаченной жидкости


где –  модуль плоской деформации, E – модуль Юнга, ν – коэффициент Пуассона, n – показать поведения жидкости, k’ – коэффициент густоты потока, Q – расход жидкости, H – мощность пласта, Ds – контраст напряжений, L – полудлина трещины, V – объем закаченной жидкости, Cl – коэффициент утечек по Картеру, K – коэффициент трещиностойкости. 

Сравнив полученные безразмерные переменные с приведенными в статье (José I. Adachi et al., 2010) для модели Pseudo3D, можно увидеть, что они идентичны с точностью до констант. Единственным важным отличием является параметр g. Из работы видно, что если поделить мощность пласта на масштабный фактор длины из этой работы, то с точностью до константы получается безразмерный параметр g. Его физический смысл – отношение мощности пласта к достигаемым длинам трещины. 

В итоге, обезразмеривание понижает размерность задачи на 5, так как модуль Юнга, коэффициент Пуассона, контраст напряжений, мощность пласта и вязкость жидкости заменяются одним безразмерным параметром g (7). 

Получение безразмерных параметров позволяет провести анализ уже проведенных работ ГРП, а именно, выявить некоторые закономерности при сравнении полудлин, полученных по дизайну и по ГДИС, в сопоставлении с общим объёмом закачки. Таким образом, необходимыми данными для анализа помимо значений полудлин трещин являются: мощность пласта, модуль Юнга и коэффициент Пуассона (входят в модуль плоской деформации), контраст напряжений в пласте, реология закачиваемой жидкости, расход жидкости и общий объём закачки жидкости ГРП. Прочностные свойства горных пород в данном анализе не учитываются. На рис. 5-7 представлены зависимости безразмерной длины от безразмерного объема для трех типов закачиваемой жидкости. 

Из графиков видно, что прогнозируемая степенная зависимость безразмерных параметров сохраняется, различается лишь степень, зависящая от реологии закачиваемой жидкости. При этом степенная зависимость сохраняется как для данных, полученных по дизайну, так и ГДИС. Приведенные выше результаты доказывают, что зависимость безразмерной длины трещины от безразмерного объема – степенная:


тогда подставив выражения (7) и (8) в (12), можно получить формулу, отражающую зависимость размерной длины трещины от объема закачки для произвольной степени a: 

Из графиков (Рис. 5-7) можно заключить, что степень зависимости безразмерной длины трещины от объема закачки a находится в диапазоне 0.6-0.7. Приняв степень зависимости равной 0.6, а показатель поведения жидкости n равным 0.5, можно получить эмпирическую формулу для расчета полудлины трещины: 

Данная формула позволяет оценить влияние ошибки (геомеханические параметры) или изменения входных параметров (технологические параметры) на изменение значения полудлины трещины (Табл. 1).

Табл. 1. Влияние изменения параметров на изменение длины трещины

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


Рис. 5. Зависимость безразмерной длины трещины от безразмерного объема: жидкость 1, м/р A; α = 0.69


Рис. 6. Зависимость безразмерной длины трещины от безразмерного объема: жидкость 2, м/р A; α = 0.63


Рис. 7. Зависимость безразмерной длины трещины от безразмерного объема: жидкость 3, м/р A; α = 0.66

Литература 

  1. Байков В.А., Бакиров Н.К., Яковлев А.А. (2010). Новые подходы в теории геостатистического моделирования. Вестник УГАТУ, 37(2), с. 209-215. Байков В.А., Бакиров Н.К., Яковлев А.А. (2012). Математическая геология. Том 1. Введение в геостатистику. М.-Ижевск: ИКИ, Библиотека нефтяного инжиниринга, 228 с. 
  2. Дьячук А.И. (2003). Некоторые аспекты реконструкции систем сбора продукции скважин на поздних стадиях разработки нефтяных месторождений. Проблемы сбора, подготовки и транспорта нефти и нефтепродуктов, 62, с. 87-97. 
  3. Дюбрюль О. (2002). Использование геостатистики для включения в геологическую модель данных. EAGE: SEG. 295 c. 
  4. Пригарин С.М., Михайлов Г.А. (2005). Методы численного моделирования случайных процессов и полей. Новосибирск: ИВМИМГ СО РАН, 259 c. 
  5. Ситников А.Н., Пустовских А.А., Рощектаев А.П., Анджукаев Ц.В. (2015). Метод определения оптимального времени отработки нагнетательных скважин в системе разработки месторождения. Нефтяное хозяйство, 3, c. 84-87. 
  6. Хасанов М.М., Белозеров Б.В., Бочков А.С., Ушмаев О.С., Фукс О.М. (2014). Применение спектральной теории для анализа и моделирования фильтрационно-емкостных свойств пласта. Нефтяное хозяйство, 12, c. 60-64. 
  7. Хасанов М.М., Белозеров Б.В., Бочков А.С., Фукс О.М., Тенгелиди Д.И. (2015). Автоматизация литолого-фациального анализа на основе спектральной теории. Нефтяное хозяйство, 12, c. 48-51. 
  8. Хасанов М.М., Мельчаева О.Ю., Ситников А.Н., Рощектаев А.П. (2013). Динамика добычи из скважин с гидроразрывом пласта в экономически оптимальных системах разработки. Нефтяное хозяйство, 12, c. 36-39. 
  9. Хасанов М.М., Падерин Г.В., Шель Е.В., Яковлев А.А., Пустовских А.А. (2017). Подходы к моделированию гидроразрыва пласта и направления их развития. Нефтяное хозяйство, 12, c. 37-41. 
  10. José I. Adachi, Emmanuel Detournay, Anthony P. Peirce. (2010). Analysis of the classical pseudo-3D model for hydraulic fracture with equilibrium height growth across stress barriers. International Journal of Rock Mechanics and Mining Sciences, 47(4), pp. 625-639. 

Сведения об авторах 

М.М. Хасанов – доктор тех. наук, директор дирекции по технологиям «Газпромнефти», генеральный директор ООО «Газпромнефть НТЦ» Россия, 190000, Санкт-Петербург, наб. реки Мойки, д. 75-79, литер Д 

А.Н. Ситников – заместитель генерального директора по научному инжинирингу, ООО «Газпромнефть НТЦ» Россия, 190000, Санкт-Петербург, наб. реки Мойки, д. 75-79, литер Д 

А.А. Пустовских – канд. физ.-мат. наук., начальник департамента научно-методического сопровождения разработки и профессионального развития ООО «Газпромнефть НТЦ» Россия, 190000, Санкт-Петербург, наб. реки Мойки, д. 75-79, литер Д E-mail: Pustovskikh.AA@gazpromneft-ntc.ru 

А.П. Рощектаев – канд. физ.-мат. наук, ведущий эксперт департамента научно-методического сопровождения разработки и профессионального развития ООО «Газпромнефть НТЦ» Россия, 190000, Санкт-Петербург, наб. реки Мойки, д. 75-79, литер Д 

Н.С. Исмагилов – канд. физ.-мат. наук, руководитель направления Управления комплексного проектирования разработки, ООО «Газпромнефть НТЦ» Россия, 190000, Санкт-Петербург, наб. реки Мойки, д. 75-79, литер Д

Возврат к списку