Одной из тенденций последних лет является повышение роли геомеханического моделирования при разработке месторождений [1]. Действительно, геомеханические процессы зачастую вносят существенные корректировки в программы освоения месторождений. К таким процессам, помимо классических задач по обеспечению устойчивости ствола скважины при бурении, относятся явления формирования трещин гидроразрыва пласта [2], уплотнения коллекторов [3, 4], реактивации крупномасштабных разломов и микромасштабных систем трещиноватости [5], эффекты самопроизвольного роста трещин гидроразрыва пласта (автоГРП) [6] и многие другие. Некоторые явления могут быть косвенно учтены посредством введения дополнительных, как правило, корреляционных зависимостей в гидродинамическом моделировании. Например, использование табулированных множителей на пористость и проницаемость в зависимости от пластового давления (RockTab) позволяет косвенно описать эффекты уплотнения коллектора в гидродинамической модели. Однако ряд явлений таким способом описать крайне сложно, например, прогнозирование эффекта автоГРП. К тому же при попытке применить корреляционный подход на других месторождениях или даже участках месторождений исследователь столкнется с необходимостью корректировки используемых корреляций. Так, хорошо известно об использовании корреляций для прогнозирования давления смыкания с использованием данных о пластовом давлении [7]. Здесь корреляции строятся отдельно на каждый продуктивный пласт. При этом построение корреляций требует каждый раз новых исследований.
Таким образом, более целесообразным представляется подход, при котором все возникающие задачи разработки и бурения решаются на базе единой модели механических свойств и напряженно-деформированного состояния. Для построения таких моделей в условиях сложной геометрии месторождения или изменчивых механических свойств необходимо использование полноценного численного расчета полей напряжения и деформации в трехмерной (3D) постановке.
В рамках политики импортозамещения в 2022 году ПАО «НК «Роснефть» начала проект по разработке собственного 3D/4D геомеханического симулятора «РН-СИГМА». В настоящее время закончена первая стадия проекта: разработано вычислительное ядро — интерфейс для его полноценного использования в ПК «РН-СИГМА». На текущий момент завершается этап апробации и тестирования 3D-ядра. В настоящей работе приведены первые результаты расчетов для участка пласта ПК1 верхнепокурской подсвиты Северо-Комсомольского месторождения, использованные для практических целей.

Геомеханическое моделирование в зависимости от постановки задачи и геологических условий можно разделить на три класса.
Одномерное геомеханическое моделирование (1D ГМ) по большей части основано на аналитических решениях задачи о напряженно-деформированном состоянии пласта или ствола скважины. 1D ГМ подразумевает расчет механических свойств пород в виде профилей вдоль ствола скважины. Класс решаемых задач включает моделирование устойчивости ствола скважины и элементов ее конструкции во время бурения или в процессе эксплуатации, устойчивость разломов и систем трещиноватости в пластах простой геометрии, подготовку геомеханических моделей для задач моделирования гидроразрыва, определение условий начала выноса песка. 1D ГМ применяется в случае слаборасчлененных пластов с границами, почти параллельными дневной поверхности. Практическое применение 1D ГМ — проектирование и сопровождение бурения скважин, проектирование дизайнов гидроразрыва пласта, оценка предельных условий безопасной эксплуатации скважин.
Трехмерное геомеханическое моделирование (3D ГМ) решает тот же класс задач, но в условиях произвольной расчлененности и геометрии залегания, в том числе при наличии разного рода тектонических нарушений. Среди дополнительных задач можно выделить задачи оценки степени уплотнения коллектора и проседания дневной поверхности. Модели расчета напряженно-деформированного состояния можно разделить на упругие и упругопластические. Последние позволяют также провести оценку возможности необратимого пластического уплотнения коллекторов или дилатантного разуплотнения с образованием зон кинжальных прорывов. Отдельным классом можно выделить задачи устойчивости сопряжения стволов скважины (fishbone).
Четырехмерное геомеханическое моделирование (4D ГМ) решает сопряженные задачи гидродинамики и геомеханики в масштабе месторождения. 4D ГМ позволяет рассчитывать изменение во времени свойств модели и состояние отдельных ее элементов, например, для учета влияния на общее напряжено-деформированное состояние истории изменения порового давления или водонасыщенности. При этом выделяют однонаправленное и двунаправленное сопряжения. Однонаправленным называют сопряжение, при котором учитывается только воздействие гидродинамических процессов на механическое состояние, но не наоборот. По сути, этот способ моделирования можно отнести к последовательному применению 3D ГМ. Двунаправленным называют сопряжение, при котором механические и гидродинамические процессы влияют друг на друга. В основном 4D ГМ применяется для полноценного моделирования гидродинамики с учетом изменения механического состояния горных пород, решения динамических задач устойчивости коллекторов, содержащих системы трещиноватости, а также для прогнозирования изменения коллекторских свойств в процессе разработки.
В настоящей работе рассматривается задача разработки ядра для 3D ГМ в упругой постановке.
В дальнейшем планируется доработка ядра для максимального охвата геомеханического моделирования.

Процесс построения трехмерной (3D) геомеханической модели можно условно разделить на два блока: построение модели свойств и расчет напряженно-деформированного состояния. Первый блок выполняется на геологической сетке с применением стандартных методов распространения свойств в межскважинном пространстве и по сути своей мало чем отличается от распространения, к примеру, фильтрационно-емкостных свойств. Для этой цели допустимо применение любых стандартных геологических или петрофизических алгоритмов (интерполяционных, статистических и др.), реализованных на корпоративной платформе RexLab ПАО «НК «Роснефть», в том числе в программных продуктах «РН-ГЕОСИМ» и «РН-СИГМА». В данной работе мы не будем детально останавливаться на этих алгоритмах.
Основной интерес представляет второй блок — расчет напряженно-деформированного состояния. Это принципиально новый алгоритм в сравнении с одномерным моделированием. В его основе лежит численное решение задачи упругости. Сформулируем центральную задачу.
Массив горной породы, включающий месторождение, находится на глубине в условиях действия на него сил со стороны вмещающей среды. Под вмещающей средой, как правило, понимают достаточно большой участок геологической формации, такой, что от его размеров при заданных внешних воздействиях напряженно-деформированное состояние интересующего массива горной породы практически не зависит. Более подробно на свойствах вмещающей среды остановимся позднее. Массив горной породы сложен из геологических сред с различными упругими свойствами. Каждую геологическую среду будем считать пористым насыщенным телом. Распределение необходимых упругих характеристик — модуля Юнга E, коэффициента Пуассона v, коэффициента
Био α, а также плотности ρ и коэффициента линейного температурного расширения β — будем считать заданными во всем массиве горной породы.
Уравнение равновесия и уравнение состояния изотропной пористой насыщенной среды с учетом термических напряжений имеют классический вид [8]:

Здесь напряжения и деформации имеют положительное значение при сжатии, в соответствии с принятым в рудной механике соглашением.
Важным вопросом при построении модели является постановка корректных граничных условий. Рассмотрим этот вопрос подробнее. На рисунке 1 изображен общий вид сеточной модели массива горной породы. Здесь можно выделить три различных типа границы: ω_up — верхняя грань, ω_side — набор боковых граней и ω_down — нижняя грань.

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

где t — вектор сил, действующих на границу; n — нормаль к границе; σ_V — вес горной породы, лежащей над верхней гранью ω_up или горное давление.

Боковые грани ω_side также являются гранями с заданными нагрузками. Граничное условие на боковых гранях записывается в форме:

Здесь σ_ij^reg называется тензором регионального напряжения, который определяется физическими процессами, происходящими на масштабах, гораздо больших масштаба выделенного массива горных пород. Как правило, эти процессы обусловлены движением геологических плит. Общим местом во всех подходах является способ определения вертикального напряжения как интегрального веса горной породы на глубине z. Касательные напряжения в плоскости x0y принимаются равными нулю:

Силы, действующие на боковые грани, далее считаем действующими в плоскости x0y. Соответствующий этим силам тензор напряжения можно характеризовать двумя главными значениями напряжения σ_h,H^reg, которые называются главными горизонтальными напряжениями, а также азимутом φ_H, определяющим ориентацию максимального главного напряжения в пространстве относительно направления на север.
При определении главных горизонтальных напряжений используются три наиболее распространенные модели.
1. Пропорциональное нагружение — способ, при котором горизонтальные напряжения пропорциональны горному давлению:

где k_h,H — заданные константы.

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

где σ_h,H,0 — горизонтальные напряжения на определенной глубине z = z₀; z=z₀; h_h,H — градиенты горизонтальных напряжений.

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

где ε_h,H,0 — тектонические деформации на определенной глубине z = z₀; ϑ_h,H — градиенты тектонических деформаций.

Кроме главных напряжений также задается азимут действия максимального горизонтального напряжения φ_H, после чего тензор регионального напряжения приводится к системе координат, в которой определена сетка, путем простого поворота вокруг вертикальной оси.
Наконец нижняя грань ω_down представляет собой абсолютно жесткий фундамент, деформация которого невозможна, по крайней мере, в нормальном направлении. Однако и тут можно рассмотреть два типа граничных условий.
1. Полное закрепление подразумевает отсутствие смещения фундамента в любом направлении:

2. Нормальное закрепление подразумевает отсутствие смещения фундамента в нормальном направлении и отсутствие касательных сил в продольном направлении:

Процесс моделирования включает несколько стадий. Первая стадия — наращивание геологической сетки. Этот процесс необходим в первую очередь потому, что описанные граничные условия применимы на масштабах, гораздо больших масштабов моделируемого массива горной породы. Кроме того, конечный размер модели приводит к тому, что узлы и ребра становятся искусственными концентраторами напряжения, что приводит к искажению решения на достаточно большом расстоянии от этих элементов. Наращивание модели снижает влияние этих паразитных напряжений и способствует повышению точности постановки граничных условий. Далее всюду будем различать исходную или геологическую сетку и наращенную сетку (геологическую сетку с наращенной вмещающей средой).
Следующая стадия представляет собой инициализацию физических свойств. На этой стадии возможно как применение уже рассчитанных свойств в исходной области, так и пересчет модели на новой сетке с использованием одномерных скважинных моделей механических свойств и стандартных алгоритмов распространения свойств. Вмещающая среда, как правило, характеризуется отсутствием данных об упругих свойствах. Распространение свойств во вмещающей среде возможно с применением алгоритмов распространения с учетом трендов. Однако при отсутствии и такой информации возможно определение упругих свойств как средних свойств исходного массива горной породы. Такой способ позволяет избежать возникновения существенных скачков напряжений и деформаций непосредственно на границе вмещающей среды и исходной модели. Для задания областей определения материалов с постоянными свойствами в ПК «РН-СИГМА» реализован отдельный механизм.

Единственным параметром, к определению которого стоит подойти более скрупулезно, является плотность горной породы над областью исходной модели. Именно плотность будет определять величину горного давления — основного источника напряженного состояния массива горной породы. Для определения этого параметра стоит воспользоваться восстановленными в рамках построения одномерных моделей трендами плотности до поверхности.
После инициализации свойств на наращенной сетке необходимо задать распределения (кубы) порового давления и температуры. Они могут быть заданы как внешние загруженные данные на исходной сетке. В простейших случаях (линейная зависимость от глубины) эти кубы могут быть рассчитаны с использованием инструментов ПК «РН-СИГМА».
Завершающий этап — выбор граничных условий из набора, описанного выше. Необходимо отметить, что на каждом этапе может быть создано произвольное количество вариантов для моделирования. Соответственно каждый расчет представляет собой конструктор, собранный из различных вариантов инициализации физических свойств, внешних полей давления и температуры и граничных условий. Это позволяет рассчитывать несколько вариантов моделей на заданной сетке и в дальнейшем выбирать наиболее подходящий для решения поставленной задачи.

Для решения поставленной задачи было реализовано вычислительное ядро на основании алгоритма конечно-элементного типа с использованием элементов первого порядка [8]. Вычислительное ядро было также внедрено в оболочку ПК «РН-СИГМА». Ядро работает с основным типом сеточной модели, используемой в геологическом и гидродинамическом моделировании — геометрией угловой точки. Единственное существенное требование к сеткам при конечно-элементном моделировании — это отсутствие разрывов. В случае, если исходная сетка содержит разрывы, они стыкуются на этапе наращивания вмещающей среды.
Сначала рассмотрим результаты сравнения расчетов с использованием разработанного ядра с результатами расчетов, выполненных в ПО «ANSYS». Геометрия задачи — однородный куб со стороной 7 м со сферическим включением диаметром 1 м. Куб находится под действием собственного веса, дополнительные боковые силы отсутствуют, смещение нижней грани равно нулю во всех направлениях. Свойства вмещающей среды: E = 7 ГПа, v = 0,45, ρ = 2 700 кг/м³; свойства включения: E = 70 ГПа, v = 0,1, ρ = 2 700 кг/м³, эффекты пороупругости отсутствуют. Свойства сред намеренно выбраны исключительно контрастные. Рассчитанное поле смещения и относительная разница между результатами расчетов представлена на рисунке 2. Видно, что различие между расчетами составляет не более 0,45 % во всей области.

Далее рассмотрим аналитический тест: решение задачи упругости в модели ограниченной деформации. Геометрически (рис. 3) модель представляет собой включение в форме параллелепипеда во вмещающей среде. Сечение включения в плоскости x0y — квадрат со стороной 1 000 м, размер по вертикали — 20 м, верхняя грань расположена на глубине 400 м от поверхности. Свойства вмещающей среды: E = 65 ГПа, v = 0,2, ρ = 2 720 кг/м³, свойства кубического включения:
E = 15 ГПа, v = 0,25, ρ = 2 650 кг/м³, пороупругие эффекты отсутствуют. Для нагружения среды использована модель заданной тектонической деформации с параметрами: ε_H = 5×10^-4, ε_h = 2×10^-4, φ_H = -15°. В такой постановке напряженное состояние вдоль линии, проходящей через центр, описывается аналитическим решением (7).

Результаты сравнения компонент главных напряжений, рассчитанных по аналитической модели, а также с использованием ПК «РН-СИГМА» и одного из промышленных аналогов (ПО «Petrel Geomechanics», «VISAGE»), приведены на рисунке 4. Можно отметить хорошую сходимость решения к результатам аналитического решения, ошибка составляет порядка 5 %.
Кроме того, видно, что отличия результатов расчета от аналитического решения для обоих вариантов симуляторов находятся на одном уровне и качественно повторяют поведение друг друга. Различие результатов решения задачи между собой также находится на достаточно низком уровне ~2–3 % в сравнении с аналитическим решением.

Применение 3D геомеханического моделирования особенно актуально в двух случаях: формации с высококонтрастными свойствами по латерали и существенное отличие геометрии пластов от плоскопараллельного чередования слоев, например, вызванное наличием разломов. Возможность выполнения расчетов на таких моделях представляет особенный интерес. Поэтому для верификации на реальных данных был выбран один из участков пласта ПК-1 Верхнепокурской подсвиты Северо-Комсомольского месторождения. Данные отложения представлены рыхлыми песками и слабосцементированными песчаниками с многочисленными тонкослоистыми и часто прерывистыми прослоями глинистых пород. Отдельные пропластки характеризуются наличием включений с высоким контрастом механических свойств. Кроме того, выбранный участок пересекает несколько разломов.
При построении 3D геомеханической модели целевых отложений был выполнен ряд шагов: проанализированы результаты керновых исследований, данные геофизических исследований (ГИС), включая специальный комплекс, результаты исследований на предварительной стадии операции гидроразрыва пласта (мини-ГРП), а также данные бурения. В результате выбраны три опорные скважины с максимальным набором данных для построения 1D геомеханических моделей. Калибровка исходных 1D геомеханических моделей выполнена на данные мини-ГРП и зафиксированные буровые события. Оценка порового давления в опорных скважинах получена на основе интерпретации данных гидродинамического каротажа. В качестве модели напряженного состояния для 1D-расчета была выбрана модель с тектонической деформацией. Значения главных тектонических деформаций оценены посредством калибровки профиля минимального горизонтального напряжения на данные о давлении смыкания по мини-ГРП и составили:
ε_h = 1×10^-4, ε_H = 2×10^-4. Коэффициент Био α по результатам адаптации принят равным 0,9 в интервалах песчаника и 1 в интервалах глин. В целом соответствие построенных 1D геомеханических моделей с имеющимися калибровочными данными оценено как высокое.

На основе действующей геологической модели в пределах зон целевых кустовых площадок построены секторные 3D геомеханические модели механических свойств согласно описанной выше последовательности действий. Сначала выполнено наращивание модели вверх до дневной поверхности, вниз до заданной глубины (1 200 м) и в горизонтальной плоскости по 200 м с каждой стороны. Далее для исключения негативного влияния сетки на результаты численного моделирования проведен тест на однородном материале. После этого выполнено распространение механических свойств пород на основе интерполяционных алгоритмов и с учетом действующей литолого-фациальной модели целевых отложений. В качестве вмещающей среды выбрана среда со свойствами, близкими к средним свойствам целевого объекта. Далее рассчитан куб начального пластового давления. На финальном шаге выполнен расчет напряженно-деформированного состояния пласта и проведена калибровка 3D геомеханической модели на результаты 1D геомеханического моделирования опорных скважин. На рисунке 5 приведен пример расчета модели выбранного участка (куб минимального горизонтального напряжения).

На рисунке 6 приведен результат верификации главных напряжений для моделируемого сектора: сравнение результатов 3D- и 1D-моделирования для одной из опорных скважин. По результатам расчетов калибровка прошла успешно: данные 1D и 3D хорошо соотносятся между собой, различие в расчетах не превышает 10 %. Исключение составляет плотный прослой в нижней части коллектора. Это связано с двумя причинами. Во-первых, данный прослой в модели не представлен отдельным слоем, следовательно, его механические свойства смешаны с выше- и нижележащими интервалами. Во-вторых, анализ геологической модели показывает, что плотные прослои на данном участке не являются протяженными, следовательно, одномерные модели не могут адекватно описывать напряжения в них.

Общий размер модели составляет примерно 1,8 млн ячеек (137×106×126), размер ячеек —
100×100 м. Время расчета модели — порядка 18 минут на современном персональном
компьютере (процессор Intel i7-8700 с частотой 3,2 ГГц, оперативная память 32 Гб).
Результаты выполненной работы по расчету 3D ГМ в настоящее время используются для геомеханического сопровождения бурения отложений пласта ПК-1, а также перекрывающих глинистых отложений Березовской и Кузнецовской свит. На основе 3D-модели построены предбуровые модели устойчивости ствола скважин: выданы рекомендации по статической плотности бурового раствора и предельно допустимой эквивалентной циркуляционной плотности раствора, определены основные зоны нестабильности ствола скважин и зоны возможных поглощений раствора. Всего за период 2022 года предбуровые модели построены для 44 скважин. По результатам бурения выполняется постбуровой анализ: на основе фактических данных ГИС имеющиеся предбуровые модели проверяются на корректность. Выполненные расчеты устойчивости ствола на основе данных 3D и на основе фактических данных ГИС, записанных во время бурения, показывают, что в целом по разрезу обе модели дают удовлетворительную сходимость с фактическими данными кавернометрии за исключением локальных пропластков. Также результаты 3D ГМ будут использованы в качестве основы для расчетов в рамках двунаправленной 4D ГМ для оценки влияния геомеханических эффектов на фильтрационные процессы в пласте при эксплуатации скважин.