CC BY
78
МАТЕМАТИЧЕСКИЕ МОДЕЛИ ДЛЯ ФИЛЬТРАЦИОННО-ПРОЧНОСТНОГО РАСЧЕТА ПРИЗАБОЙНЫХ ЗОН СКВАЖИН
А.Л. Ковалёв, Е.В. Шеберстов (ООО «Газпром ВНИИГАЗ»)
Выбор конструкции забоя, обеспечивающей высокую производительность скважины и исключающей разрушение коллектора, представляет собой задачу, для решения которой на месторождениях с большим эксплуатационным фондом можно использовать эмпирический подход. Для прогноза прочности вскрываемого разреза морских месторождений, экономичная разработка которых требует применения небольшого числа горизонтальных скважин, при выборе способа их заканчивания приходится максимально использовать данные геофизики и испытаний керна. Примером может служить подход к исследованию риска выноса песка на скважинах Штокмановского ГКМ [1]. Для выполнения такой программы необходимо специальное испытательное оборудование и математические (в том числе компьютерные) модели, позволяющие решать следующие задачи:
■ подбор и адаптация математических моделей деформирования горных пород по результатам стендовых испытаний;
■ построение фильтрационно-прочностной модели для прогноза напряжений на различных участках месторождения;
■ построение фильтрационно-прочностных моделей прискважинных участков пласта и прогноз их состояния.
Методической основой рассматриваемого подхода служит предположение о том, что заключенный в керне материал сохраняет свойства, которыми определяется поведение горной породы в окрестности ствола скважины и что выявить эти свойства можно, выполнив испытания образцов по специальной программе. Обычно модель материала выбирают на основании трехосных испытаний, а ее адекватность проверяют по испытаниям в условиях, приближенных к реальным.
Трехосные испытания проводят на образцах кубической или цилиндрической формы. Кубические образцы (сторона 40 или 50 мм) из керна Штокмановского и Приразломного месторождений испытывались в Институте проблем механики РАН на испытательном стенде трехосных независимых нагружений (ИСТНН) [1]. Согласно общепринятой схеме, на первом этапе создается гидростатическая (одинаковая для всех граней) на-
грузка, постепенно увеличивающаяся до определенного значения (давление обжима). Затем нагрузка на боковые грани поддерживается постоянной, а нагрузку на горизонтальную грань продолжают увеличивать до появления признаков разрушения образца.
Визуальный анализ показывает, что на начальном участке (от атмосферного давления до давления обжима, а иногда и выше) кривые деформирования не прямолинейны. Эту особенность поведения горных пород отмечают многие исследователи. Она объясняется тем, что деформирование на начальном участке происходит вследствие закрытия микротрещин, появившихся при извлечении керна из пласта. Прямолинейный участок появляется при более высоких нагрузках и завершается участком с увеличенной скоростью деформации при постоянной или даже снижающейся нагрузке (рис. 1).
При выборе модели материала предполагалось, что образец однороден и главные напряжения в каждой точке образца совпадают с напряжениями на гранях. Это же относится и к деформациям. Начальный участок нелинейного поведения не учитывался, т.к. для пластовых условий характерны более высокие осевые напряжения. Участок прямолинейной зависимости представлялся моделью изотропного пороупругого тела (модель Гука - Терцаги - Био).
где о. - эффективные напряжения; X, О - коэффициенты Ламе; К - объемный модуль; £у - тензор деформаций (е = 8.); в - коэффициент сжимаемости зерен скелета; 5. - символ Кронекера. Поскольку испытания проводили в атмосферных условиях, то при адаптации последний член формулы (1) приравнивался нулю. При выходе напряжений на поверхность течения, определяемую уравнением
наряду с упругими (обратимыми) деформациями ее., вычисляемыми по (1), возникают пластические (остаточные) деформации ер:
8.. = £е. +£р . (3)
У У У V /
При / (о., к) < 0 материал находится в упругом состоянии (здесь к -параметр упрочнения).
Определяющие соотношения в зоне пластичности заданы в виде
Сту =Х е-5. +20 • 8у +р5 • К - РЬ у,
(1)
/ (Су, к) = 0,
(2)
где V - тензор скорости деформации; gо) - потенциал течения. Если g = / то (4) представляет так называемый ассоциированный закон течения. В программах, созданных в Центре исследований нефтегазовых пластовых систем и технологического моделирования, в качестве поверхности течения использованы критерии Кулона - Мора (/=/с - м) и Друкера - Прагера
где о1 и о3 - максимальное и минимальное главные эффективные напряжения соответственно; 11 и 32 - первый инвариант тензора напряжений и второй инвариант девиатора тензора напряжений. Константы к, С, а, К выражаются через когезию (8С) и угол внутреннего трения (ф).
Выбор моделей производился в следующей последовательности. Для каждого интервала разреза проводилось три испытания с давлением обжима 2, 10, 20 МПа. По трем кругам Мора известным образом определялись когезия и угол внутреннего трения, затем вычислялись коэффициенты в моделях (5, 6). В процессе испытания производилась цифровая запись нагрузок и деформаций. С помощью выбранной модели материала для фактической последовательности напряжений определялись деформации, и расчетные значения сопоставлялись со значениями, измеренными в опыте. Вычисление деформаций производилось по (1), если образец находился в упругом состоянии. В противном случае строилась упруго-пластическая матрица, связывающая приращения деформаций и напряжения, и определялись пластическая и упругая составляющие деформаций. Были рассмотрены две модели: идеально-пластическое деформирование / (о., к) = =/(оу) и модель с упрочнением. Поверхность течения описывалась критерием Друкера - Прагера. Параметр упрочнения определялся накопленной пластической деформацией вдоль оси 2
где 8^ - когезия в состоянии перехода от упругих к пластическим дефор-
кущая пластическая деформация по оси 2 (как оси наибольшего напряжения в трехосных испытаниях); параметры 8е и ер2тах подбирались в про-
/а -м - °1 - Ь 'а з - С ; їв-м - а ' 11 + 4^2 - К»
(5)
(6)
через которую выражалась когезия:
2 тах
(7)
мациям; 8реак - когезия в состоянии максимального напряжения; ер - те-
цессе адаптации. Второй параметр использовался в прогнозных расчетах в качестве критерия прочности. В принятой модели угол внутреннего трения постоянен. Пример воспроизведения моделью второго цикла трехосного испытания (давление обжима 10 МПа) образца из разреза Приразломного месторождения (карбонатный коллектор) приведен на рис. 1.
Добиться хорошего соответствия расчетных и опытных данных для образцов из керна Приразломного месторождения удавалось не всегда. Как правило, проблемы возникали, когда три круга Мора плохо описывались прямолинейной огибающей. Можно предположить две причины этих расхождений: анизотропия породы и неоднородность материала даже в пределах керна, из которого изготавливались «одинаковые» образцы.
Для проверки адекватности моделей использованы результаты двух экспериментов (выполненных также на установке ИСТНН). В первом моделировались условия на стенке скважины при вскрытии пласта. Модель идеальной пластичности, адаптированная по трехосным испытаниям, хорошо воспроизвела переход в состояние разрушения.
Второй эксперимент проводился по схеме «полый цилиндр». Испытывались образцы кубической формы с отверстием. Создавалась гидростатическая нагрузка на внешние грани (включая торец с отверстием), а стенка отверстия контактировала с атмосферой. Начиная с некоторой нагрузки, с поверхности отверстия отделялись частицы, накопленный вес которых определялся на весах, расположенных непосредственно под отверстием. В этих экспериментах прочность образцов оказалась существенно выше модельной. Подобный эффект хорошо известен из литературы, но его механизм до сих пор не объяснен. Одна из возможных причин -влияние масштабного фактора.
Модели материала, полученные на основании трехосных испытаний, используют для создания фильтрационно-прочностных моделей участков пластовой системы. Поскольку давление флюида и напряженно-деформированное состояние (НДС) пласта в межскважинном пространстве слабо зависит от деталей конструкции забоя, фильтрационно-прочностной расчет можно проводить в два этапа. На первом этапе используется полномасштабная модель, с помощью которой производится «глобальный» расчет разработки месторождения на несколько десятилетий с использо-
о о г\ _
ванием одной из возможных конструкций скважин. Затем, используя результаты глобального расчета, можно корректировать конструкцию забоя с помощью локальных моделей.
Образец П-5 цикл 2
-8-4 -2 0 2 4 в 8 10 12 14
___________________Деформация 10ОО__________________
Рис. 1. Воспроизведение второго цикла трехосных испытаний для образца П-5: ei, e2, e3 — фактические кривые по трем осям; e1m, e2m, e3m — модельные кривые (напряжение и деформация на сжатие — положительные)
Для условий Приразломного месторождения полномасштабная фильтрационно-прочностная модель была создана с использованием средств симулятора Eclipse версии 2005 г. Модель была опробована на фрагменте залежи, границы которого из соображений симметрии можно считать непроницаемыми. Для геомеханической части модели на верхней границе задано сжимающее напряжение, соответствующее давлению вышележащего массива пород, на прочих границах заданы нулевые перемещения в направлениях, перпендикулярных к границам. По мере накопления опыта моделирования геомеханических процессов вопрос о граничных условиях для фрагментов будет уточняться.
Следует отметить, что при этом типе граничных условий соотношение между величинами вертикальных и горизонтальных эффективных напряжений определяется исключительно значением коэффициента Пуассона, варьируя которым, можно добиться требуемой величины бокового распора.
Модель фрагмента пласта представлена на рис. 2. Размеры фрагмента по осям: X - 2000 м; Y - 1000 м; Z - 50 м; размерность сетки - 37 х 10 х 11 блоков. Фрагмент пласта дренируется двумя половинами горизонтальных скважин по краям и одной полной горизонтальной скважиной в центре.
Горизонтальные стволы скважин направлены вдоль оси У и полностью вскрывают фрагмент на этом направлении. Размеры блоков сетки, содержащих горизонтальные стволы, в перпендикулярной им плоскости 2 х 2 м. В целом такая схема соответствует усредненной схеме проектного размещения скважин. Изначально фрагмент полностью насыщен нефтью, что отражает состояние в центральной сводовой части залежи, где имеется сухое поле.
Расчеты с использованием указанной модели были выполнены для двух вариантов: разработка с поддержанием пластового давления (ППД) посредством закачки воды (рекомендовано проектом) и разработка в режиме истощения. При средней величине коэффициента Пуассона 0,22 (получена на керне) упругое состояние во всех сеточных блоках модели сохраняется для обоих рассмотренных режимов. При значении коэффициента 0,1 (минимальное значение) для режима разработки залежи с ППД пласт по-прежнему будет сохранять упругие свойства, тогда как при истощении пласт частично перейдет в пластическое состояние. На рис. 3 показана расчетная динамика нормальных эффективных напряжений в сеточном блоке с добывающей скважиной для варианта с ППД при коэффициенте Пуассона 0,1.
Рис. 2. Полномасштабная модель фрагмента горизонта I Приразломного месторождения
Расчеты показали, что наибольшее напряжение материала, которое определяется прежде всего разницей между эффективными нормальными напряжениями вдоль осей 2 и X, возникает в блоках с добывающей скважиной в момент наибольшего падения давления. Для разработки залежи в режиме ППД такой момент приходится на конец десятого года разработки и определяется временем прохода вала компенсации давления от забоя нагнетательной скважины к забою добывающей скважины.
Как видно из рисунка, расчетные изменения напряженного состояния в межскважинном пространстве на всем протяжении разработки с ППД невелики. Более крутые поля получаются при моделировании режима истощения.
Зная динамику напряженного состояния в межскважинном пространстве, можно перейти к анализу прочности прилегающей к скважине области пласта средствами математического моделирования. При этом приходится рассматривать упруго-пластическое поведение горной породы совместно с течением флюида. Получить аналитические решения удается только для сферической и осевой симметрий и для простых определяющих соотношений. Так, многими авторами предложены аналитические решения для притока к сферическим и цилиндрическим кавернам в рамках теории идеальной пластичности [2, 3]. Чтобы расширить область исследуемых эффектов (нелинейная многофазная фильтрация в неоднородной по радиусу среде), в Центре созданы численные модели для сферической и цилиндрической геометрий. Математические модели в данном случае сводятся к системам обыкновенных дифференциальных уравнений, для решения которых использован метод Рунге - Кутта.
Для асимметричных условий (например, приток к горизонтальной скважине при коэффициенте бокового распора < 1) необходимо прибегать к более сложным моделям. В модели Центра принята схема плоского деформированного состояния и ассоциированный закон пластического течения с упрочнением в соответствии с билинейной моделью, использованной при адаптации трехосных испытаний. Для численного решения применен метод конечных элементов, в соответствии с которым на каждом расчетном шаге выполняется минимизация функционала, выражающего принцип виртуальных перемещений [4]. Энергия деформации определяется упругой матрицей для области упругих деформаций и упругопластической - для области пластических деформаций. Внутренние силы представлены градиентом пластового давления, рассчитываемым в отдельной подпрограмме. В качестве внешних сил заданы нормальные эф-
фективные напряжения, приложенные на границах расчетной области. В области упруго-пластических деформаций задача решается в так называемой инкрементальной постановке, т.е. задаются соотношения приращения напряжений - приращения деформаций. Такая постановка позволяет учесть путь нагружения.
Рис. 3. Динамика упругих эффективных напряжений в сеточном блоке с добывающей скважиной в процессе разработки фрагмента Приразломного месторождения с ППД
(коэффициент Пуассона 0,1)
Расчетная область, представляющая собой % сечения блока со скважиной размером 2 х 2 м, была поделена неравномерной сеткой на 162 треугольных элемента со 100 узлами. На внешних границах расчетной области задавались нормальные эффективные напряжения вдоль осей X и 2 (£Х и $>2 соответственно). При этом фильтрационная нагрузка задавалась в виде узловых значений пластового давления.
Возникающая в области упруго-пластических деформаций нелинейная задача решалась явным или неявным образом. В последнем случае использовалась следующая итерационная схема:
К(ип) ип -^ + К(ип) Ьи = 0;
иг‘
п+1 „ п . £ „ п+1
= и +Ьи
(8)
где К - матрица коэффициентов уравнений (в упругой области не зависит от перемещений, в упруго-пластической - пересчитывается на каждой
итерации); ^ - вектор правых частей уравнений; и - вектор перемещений; Ьи - вектор приращений перемещений; п - номер итерации.
Сходимость процесса и момент прекращения итераций определяются по заданному значению какой-либо из норм вектора приращений перемещений (например, нормы, определяемой максимальным абсолютным значением приращения перемещения в узле). В расчетах для материала образцов керна Приразломного месторождения погрешность выполнения условия пластического течения ( = 0) не превышала 0,03 МПа (в большинстве расчетов была порядка 0,01 МПа). В условиях примера обе схемы давали близкие результаты при несколько больших затратах компьютерного времени в случае неявной схемы.
Расчеты по программе Центра были выполнены для образцов П-5 и П-9 при коэффициенте бокового распора порядка 0,2. Они показали, что во всех рассмотренных случаях пластическая зона локализуется в кольце элементов, примыкающих непосредственно к стволу скважины (ширина кольца составляет порядка 2 см). Для материала образца П-5 пластическая зона более обширна и смещена к горизонтальному радиусу. Фильтрационная нагрузка увеличивает размеры этой зоны (рис. 4). Для материала образца П-9, напротив, пластическая зона смещена к вертикальному радиусу, состоит всего из одного элемента и исчезает при задании фильтрационной нагрузки. Было установлено, что причина различий обуславливается разницей в величине угла внутреннего трения. Для образца П-9 этот угол достаточно большой - 38° против 30° для образца П-5.
Расчетные пластические деформации растяжения оказались существенно выше сжимающих, поэтому в дальнейшем желательно скорректировать программу стендовых испытаний и уточнить модель материала в области растягивающих нагрузок.
Расчетные пластические деформации растяжения оказались существенно выше сжимающих, поэтому в дальнейшем желательно скорректировать программу стендовых испытаний и уточнить модель материала в области растягивающих нагрузок.
Параллельно с разработкой собственной программы Центром был выполнен ограниченный объем расчетов в программном комплексе Л№¥8 с использованием критерия пластичности Друкера - Прагера без упрочнения. Рассчитаны пары вариантов, различающихся величиной когезии: когезия, соответствующая началу пластического течения, и когезия, соответствующая разрушению. Принят ассоциированный закон течения. Использованы две модели: модель плоского деформированного состояния (как в программе
Центра) и трехмерная однослойная модель. Внешние нагрузки для модели плоского деформированного состояния заданы такие же, как и при расчетах по программе Центра. Фильтрационная нагрузка не учитывалась. Расчетная область поделена сеткой, сгущающейся к стволу скважины, на 1936 четырехугольных элементов (44 разбиения по углу х 44 разбиения по радиусу). Расчетная схема содержит 5985 узлов, из которых 2025 размещены в углах четырехугольных элементов, а остальные - на их сторонах.
Рис. 4. Схема конечных элементов образца П-5 Приразломного месторождения (элементы, перешедшие в пластику, выделены заливкой): а - без фильтрации;
б - с фильтрацией
Трехмерная однослойная модель создана с целью воспроизведения пути нагружения материала аналогично принятому в [3]: общая деформация вдоль оси ствола сохраняет значение, соответствующее нетронутому массиву. Слой толщиной 10 см располагался перпендикулярно оси ствола. Нагружение модели производилось в три шага. На первом шаге задавались напряжения на внешних границах и напряжение на стенке скважины, равное полусумме напряжений на внешних границах; на втором шаге напряжение на стенке скважины уменьшалось наполовину; на третьем - снижалось до нуля.
На рис. 5 показаны карты распределения пластических деформаций вдоль горизонтальной и вертикальной осей, перпендикулярных стволу скважины, по результатам расчетов в ЛШУЗ однослойной задачи для материала образца П-9 с когезией начала пластики.
Рис. 5. Карта изолиний пластических деформаций прискважинной области пласта со свойствами образца П-9 Приразломного месторождения (однослойная модель, когезия начала пластики): а - вдоль горизонтальной оси; б - вдоль вертикальной оси
В таблице сведены максимальные значения растягивающих и сжимающих пластических деформаций для вариантов с когезией начала пластики (идеальная пластичность), рассчитанных в ANSYS (на плоской и однослойной моделях) и по программе Центра. Как видно, расчетные величины пластических деформаций заметно различаются. Это можно объяснить различиями в детальности сеток и типов конечных элементов. Различия в двух решениях ANSYS вызваны разными условиями нагружения вдоль оси скважины. Таким образом, использование численных методов требует тщательного контроля параметров расчетной модели и условий нагружения путем сопоставления с аналитическими решениями и результатами опытов. Существенным выводом является также необходимость уточнения поведения материала при растяжении.
В настоящей статье представлен комплекс программ, используемых Центром для обслуживания всего цикла исследований НДС продуктивного пласта. Он позволяет по данным физического эксперимента настраивать прочностные модели материала пласта и проводить фильтрационнопрочностные расчеты режимов разработки месторождения на полномасштабной и локальных моделях. Комплекс базируется как на использовании коммерческого программного обеспечения (Eclipse, ANSYS), так и на программах собственной разработки. К числу последних можно отнести программы, реализующие новые алгоритмы, в частности численные решения фильтрационно-прочностных задач для сферической и цилиндрической областей.
Максимальные значения растягивающих и сжимающих пластических деформаций (б/р) для вариантов с когезией начала пластики, рассчитанных в программном комплексе ЛМ8У8 и по программе Центра
Показатель Значение по образцам
П-5 П-9
Деформация разрушения 0,002 0,0015
Максимальная растягивающая пластическая деформация
- Программа Центра 0,000167 0,000214
- ЛNSYS (2Б) 0,001063 0,001398
- ЛNSYS (3Б) 0,003242 0,005395
Максимальная сжимающая пластическая деформация
- Программа Центра 0,000060 0,000011
- ЛNSYS (2Б) 0,000425 0,000483
- ЛNSYS (3Б) 0,001133 0,000718
Вместе с тем представленный комплекс, несомненно, нуждается в совершенствовании и расширении имеющихся функциональных возможностей. Необходимо повысить точность фильтрационно-
прочностных моделей как в части их адаптации к данным физического эксперимента, так и в части прогноза НДС. В частности, программный комплекс ЛШУ8 неудобен для решения рассмотренных задач из-за трудностей учета фильтрационных нагрузок.
Необходимо увеличить набор прочностных моделей материала пласта, увязав его с постановкой дополнительных физических экспериментов, что позволило бы учитывать в прогнозных расчетах такие факторы, как дилатансия, прочность материала на растяжение и т.п.
Необходимо расширить функциональность моделей, что дало бы возможность учитывать взаимодействие хвостовика с породой пласта; влияние обводнения, солянокислотных обработок, гидроразрыва и прочих воздействий на призабойную зону; прогнозировать образование каверн и вынос породы из скважины и т. п.
Необходимо, наконец, в рамках полномасштабного фильтрационно-прочностного моделирования пластовых систем - месторождений и ПХГ - решать сопряженные геодинамические задачи подвижек земной поверхности, сохранения целостности покрышки при нагнетании в пласт флюидов и т.п.
Список литературы
1. Климов Д.М. Определение прочностных характеристик пород Штокмановского ГКМ и оценка рисков выноса песка при его разработке / Д.М. Климов, Р.М. Тер-Саркисов, С.Е. Чигай и др. // Газовая промышленность. - 2010. - № 11. - С. 57-60.
2. Brath R.K. Stability and Failure of Sand Arches / R.K. Brath, R. Risnes // SPE J. - 1981. - April. - Р. 236-248.
3. Risnes R. Sand Stresses Around a Wellbore / R. Risnes, R.K. Bartli, P. Horsrud // SPE J. - 1982. - P. 883-898.
4. Зенкевич О. П. Метод конечных элементов в технике: пер. с англ. / О.П. Зенкевич. - М.: Мир, 1975. - 544 с.
