МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
имени М.В. ЛОМОНОСОВА
МЕХАНИКО-МАТЕМАТИЧЕСКИЙ ФАКУЛЬТЕТ
Кафедра газовой и волновой динамики
На правах рукописи
УДК 539.3
Захаров Павел Петрович
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ
ДЕФОРМИРОВАНИЯ И РАЗРУШЕНИЯ
НЕФТЕНОСНОГО ПЛАСТА
Специальность 01.02.04 – механика деформируемого твердого телаАвтореферат диссертации
на соискание ученой степени кандидата физико-математических наукМосква - 2011
Работа выполнена на кафедре газовой и волновой динамики механикоматематического факультета Московского государственного университета имени М.В. Ломоносова.
Научный руководитель: Доктор физико-математических наук, профессор Киселев Алексей Борисович
Официальные оппоненты: Доктор физико-математических наук, старший научный сотрудник Пшеничнов Сергей Геннадьевич Доктор физико-математических наук, старший научный сотрудник Рыбакин Борис Петрович
Ведущая организация: РГУ нефти и газа имени И.М. Губкина
Защита состоится «13» мая 2011 г. в 16 часов на заседании диссертационного совета Д 501.001.91 при Московском государственном университете имени М.В. Ломоносова по адресу: 119992, Москва, Ленинские горы, Главное здание МГУ, Механико-математический факультет, аудитория 16-10.
С диссертацией можно ознакомиться в библиотеке механикоматематического факультета МГУ (Главное здание, 14 этаж) Автореферат разослан «13» апреля 2011 г.
Ученый секретарь диссертационного совета Д 501.001.91, доктор физико-математических наук, профессор С.В. Шешенин 2
Общая характеристика работы
Актуальность темы. В инженерной практике широко используются создаваемые бурением искусственные горные выработки (полости) кругового сечения различного диаметра от нескольких сантиметров (шпуры) до метров (скважины, шахтные стволы и т.п.). При бурении большое значение имеет создание условий, при которых обеспечивается устойчивость от разрушения породы со стороны внутренней поверхности выработки. Характер разрушения горной породы во многом связан с наличием в ней структурных неоднородностей различных масштабов (пор, трещин и т.п.). В данной работе рассматривается численное моделирование динамики деформирования и разрушения горного пласта вблизи скважины при резком снятии внутрискважинного давления. Актуальность темы исследования связана с практической необходимостью прогнозировать поведение горных пластов в процессе их бурения и эксплуатации, в частности, при разработке месторождений в нефтегазодобывающей промышленности.
Цель диссертационной работы. Целью данной работы является разработка алгоритмов численного моделирования процессов необратимого динамического деформирования и разрушения горного пласта, определение зон и характера его разрушения вблизи скважины.
Научная новизна. Впервые получено численное решение в двухмерной плоской постановке задачи динамики необратимого деформирования и разрушения горного нефтегазонасыщенного пласта с использованием связанной модели повреждаемой термоупругопластической среды с двумя скалярными параметрами поврежденности. При этом определялись зоны разрушения в пласте, и проводилось явное построение макротрещин в процессе численного расчета с использованием разработанного автором диссертации оригинального вычислительного алгоритма.
Научная и практическая значимость. Созданные автором диссертации модели, оригинальные алгоритмы и программа численного расчета, кроме своей научной ценности, могут найти применение в горнодобывающей, нефтяной и ряде других отраслей промышленности. Достоверность результатов, полученных в диссертационной работе, обусловлена использованием термодинамически корректных моделей сплошных сред, фундаментальных законов механики и общепризнанных численных методов.
Работа выполнена при поддержке РФФИ (гранты № 06-01-00185а и № 09-01-00144а), а также компании Schlumberger.
Апробация работы.
Основные результаты диссертации докладывались на следующих научных конференциях и семинарах:
- Конференция «Ломоносовские чтения». Секция механики. Апрель 2008, 2009 и 2010 гг., Москва, МГУ;
- VIII Всемирный конгресс по вычислительной механике, 30 июня – 4 июля 2008 г., Венеция, Италия;
- Научная конференция «Современные проблемы газовой и волновой динамики», посвященная 100-летию со дня рождения академика Х.А.
Рахматулина, 21-23 апреля 2009, Москва, МГУ;
- Конференции «Теория и практика расчета зданий, сооружений и элементов конструкций. Аналитические и численные методы», 22 ноября 2009, Москва, МГСУ;
- IV Европейская конференция по вычислительной механике, Париж, Франция, 16-21 мая, 2010;
- XII Международный семинар «Супервычисления и математическое - Международный научный симпозиум по проблемам механики деформируемых тел, посвященного 100-летию со дня рождения А.А.
Ильюшина, 20-21 января, 2011 г., Москва, МГУ;
- Научно-исследовательские семинары кафедры газовой и волновой динамики механико-математического факультета МГУ (руководители:
академик РАН Е.И. Шемякин, академик РАН Р.И. Нигматулин);
- Научно-исследовательский семинар кафедры механики композитов механико-математического факультета МГУ (руководитель профессор Б.Е.
Победря);
- Научно-исследовательский семинар кафедры нефтегазовой и подземной гидромеханики факультета разработки нефтяных и газовых месторождений РГУ нефти и газа имени И.М. Губкина (руководитель профессор В.В. Кадет);
- Внутренние конференции и семинары компании Schlumberger, Москва, 2007-2009 гг.
Публикации. По работе имеется 9 публикаций, в том числе одна статья в журнале из перечня ВАК.
Диссертация работа состоит из введения, четырех глав, заключения и списка литературы. В работе содержится 54 рисунка, 91 библиографическая ссылка. Общий объем диссертации составляет 105 страниц.
Во введении обосновывается актуальность численного моделирования процессов деформирования и разрушения геологических пластов, содержащих в своих порах и трещинах нефть и/или газ. Определяются цели, новизна проведенных исследований, приведен список публикаций по теме диссертации, а также приведен список организаций и конференций, где докладывались основные результаты работы.
В первой главе представлена термодинамически корректная математическая модель термоупруговязкопластической среды с двумя скалярными параметрами поврежденности, которая является основой для моделирования поведения нефтегазонасыщенного пласта в диссертационной работе (Киселев А.Б. Математическое моделирование динамического термовязкоупругопластической среды // Вестн. МГУ. Сер.1. Матем.Механ. – 1998. - № 6). Первыми работами, в которых были введены скалярные параметры поврежденности были работы Л.М. Качанова и Ю.Н. Работнова (Качанов Л.М. О времени разрушения в условиях ползучести // Изв. АН СССР. ОТН. - 1958. - № 8; Работнов Ю.Н. Механизм длительного разрушения // Сб. «Вопросы прочности материалов и конструкций». - М.:
Изд-во АН СССР, 1959). Тензорные меры поврежденности впервые ввел А.А.
Ильюшин (Ильюшин А.А. Об одной теории длительной прочности // Инж. Ж.
Механ. Твердого тела. – 1967. - № 3).
Во второй главе приведена математическая постановка исследуемой задачи.
Рассматривается следующая задача о динамике геологической породы.
Горный пласт, плоскость которого перпендикулярна к оси скважины, находится в состояние покоя. Напряженное состояние пласта вблизи скважины определяется естественным горным давлением и внутрискважинным давлением. В определенный момент времени внутрискважинное давление падает на заданную величину, что приводит к динамическому деформированию пласта.
Приняты следующие основные допущения:
- параметры задачи не зависят от пространственной координаты, нормальной к плоскости пласта;
- предполагается плоская деформация пласта, задача решается в плоской двумерной постановке;
- пласт моделируется как повреждаемая термоупругопластическая среда;
- процесс деформирования считается адиабатическим.
Введем цилиндрическую систему координат: ось Oz совпадает с осью скважины, а пласт находится в плоскости Or (см. рис. 1).
Рис. 1. Цилиндрическая система координат, связанная со скважиной и Уравнения массы, импульса и энергии в цилиндрической системе координат имеют следующий вид:
Здесь и далее точка над символом обозначает материальную производную по времени; - плотность; Vr, V - компоненты вектора скорости; rr,, r компоненты тензора напряжений, который раскладывается на шаровую и девиаторную части – и Srr, S, Szz, Sr; rr,, r - компоненты тензора скоростей деформации; rr,, zz, r - пластические составляющие тензора скоростей деформации; T - температура;, - скалярные параметры поврежденности среды; c - теплоемкость при постоянных напряжениях, v модуль объемного расширения;, A - ассоциированные к модели параметры среды, связывающие термические процессы с процессами накопления повреждений.
Компоненты тензора скорости полных деформаций выражаются через компоненты вектора скорости следующим образом:
и раскладываются на упругую и пластическую составляющие:
Относительно пластической деформации принимается предположение о несжимаемости:
Система определяющих уравнений связанной модели повреждаемой термоупругопластической среды, полученная с использованием термодинамических принципов механики сплошной среды, имеет следующий вид:
Здесь символом обозначена яуманновская производная от компонент девиатора тензора напряжений (S ij ) = S ij S ij jk S jk ik ; eij - девиатор тензора части, по отношению к которым среда ведет себя как неповрежденная; K0 и 0 - объемный модуль и модуль сдвига неповрежденного материала соответственно; ij - тензор вихря, компоненты которого выражаются через компоненты вектора скорости следующим образом:
Последнее соотношение из (7) – закон Мизеса-Шлейхера, связывающий предел упругости Y0 и давление в пласте (-); c1, c2- константы материала.
Система уравнений (1) - (8) замыкается кинетическими уравнениями для параметров поврежденности, :
Здесь 0 - динамическая вязкость материала; p0 - начальное давление газа, заполняющего поры (горное давление); - показатель адиабаты газа(жидкости), заполняющего поры; 0 - начальная пористость среды; S u = S ij S ij - интенсивность девиатора тензора напряжений; Su - константа материала; H(…) - единичная функция Хевисайда.
Параметр поврежденности связан с шаровой частью тензора напряжений и интерпретируется как объемная поврежденность типа микропор. Параметр связан с интенсивностью девиатора тензора напряжений и описывает разрушение по сдвиговому типу.
В качестве критерия начала макроразрушения используется энтропийный критерий предельной удельной диссипации (Киселев А.Б., Юмашев М.В. Деформирование и разрушение при ударном нагружении.
Модель повреждаемой термоупругопластической среды // ПМТФ. – 1990. который в случае адиабатического процесса имеет следующий вид:
Здесь t* – время начала разрушения; D* - предельная удельная диссипация (константа материала), которая складывается из механической диссипации за счет необратимого пластического деформирования d M = Sij&ijp и диссипации континуального разрушения d F = 2 + 2.
Начальные условия. Изначально горный пласт покоится, поэтому в начальный момент времени t = 0 Vr(r,) = V(r,) = 0.
В расчетах в качестве начального распределения плотности было принято равномерное распределение: = (r,) = 0.
В качестве начального распределения напряжений в пласте используется решение следующей статической плоской задачи.
Рассматривается бесконечное призматическое тело с бесконечным круглым цилиндрическим вырезом (плоская деформация). На рисунке 2 показана схема приложения нагрузок к телу. Напряжения 1, 2 приложенные на бесконечности, соответствуют горному давлению в пласте, а напряжения соответствует внутрискважинному давлению; 1 < 0, 2 < 0, 3 < 0. В общем случае 1 2, таким образом начальное распределение напряжений зависит от полярного угла. Данная задача имеет решение в линейноупругой постановке и в упругопластической постановке, известной как задача Галина (Галин Л.А. Плоская упруго-пластическая задача // ПММ. – 1946. - Т. 10, Вып. 3).
Рис.2. Схема приложения нагрузок для начального распределения Решение в линейноупругой постановке получается суперпозицией трех решений (Тимошенко С.П., Гудьер Дж. Теория упругости. - М.:Наука, 1975) отдельно для напряжений 1, 2 и 3 и дается следующими формулами:
Задача Галина решается для условия пластичности Треска - Сен-Венана ( xx yy ) 2 + 4 r2 = 4 K 2. Решение дается формулами (13) в пластической области а в упругой области Здесь xx, yy, xy - компоненты тензора напряжений в декартовой системе координат; z = x+iy- комплексная координата; K - константа пластичности, предел текучести при чистом сдвиге; a - радиус выреза; - коэффициент существования решения Галина существенно, чтобы круговой вырез был полностью охвачен пластической областью, т.е. лежал внутри эллипса записывается в виде неравенства, связывающего параметры 1, 2, 3 и K :
Граничные условия. В момент времени t = 0 происходит резкое падение внутрискважинного давления с величины 3 на величину p > 0.
Поэтому на стенке скважины задается давление 3 + p. Случай, когда p = |3|, соответствует условию свободной поверхности.
В третьей главе описаны наиболее важные аспекты численного метода.
Поставленная задача решается численно явным конечно-разностным методом типа Уилкинса (Уилкинс М.Л. Расчет упругопластических течений // Вычислительные методы в гидродинамике. - М.: Мир, 1967) на лагранжевой сетке, которая движется и деформируется вместе со средой.
Моделирование разрушения материала осуществляется явным образом путем расщепление расчетной сетки (Немирович-Данченко М.М. Модель гипоупругой хрупкой среды: применение к расчету деформирования и разрушения // Физическая мезомеханика -1998.- Т.1.- №2). Если вычисленное в узле сетки значение удельной диссипации энергии превысило предельное, то материал считается разрушенным, и в данном узле происходит элементарное расщепление сетки по одному из вариантов, изображенном на рис. 3.
Рис. 3. Варианты элементарных трещин для внутреннего узла Вариант расщепления определяется из анализа вектора напряжений на смежных для данного узла ребрах. Взаимодействие «элементарных» трещин с уже существующими трещинами формирует зоны разрушения большего масштаба. В результате расщепления на вновь образовавшихся граничных ребрах задаются граничные условия: условие свободной поверхности или контактные условие в зависимости от ситуации.
Используемый в данной работе алгоритм реализации граничных условий на контактных поверхностях заключается в коррекции координат и скорости граничных узлов расчетных ячеек (Высокоскоростное взаимодействие тел / Под ред. В.М. Фомина. - Новосибирск: Изд-во СО РАН, 1999). Этот алгоритм является симметричным, в том смысле, что при коррекции координат взаимодействующие границы входят симметричным образом.
В четвертой главе приведены результаты расчетов и их обсуждение.
Задача обладает двумя осями симметрии, поэтому расчеты проводились для четверти области. На границах = 0 и = /2 при интегрировании уравнений движения учитывается вклад симметричных ячеек, так чтобы выполнялось условие V = 0. Внешняя граница расчетной области бралась достаточно далекой, чтобы не учитывать отражение от внешней границы волн, идущих со стенки скважины.
Расчеты проводились в трех разных постановках. В первой постановке в численном расчете не осуществлялось явного построения зон разрушения.
Во второй постановке зоны разрушения строились явным образом, однако не учитывались контактные взаимодействия на вновь образовавшихся свободных поверхностях в результате расщепления сетки. В третьем подходе на свободных поверхностях контактные взаимодействия учитывались.
В расчетах использовались следующие значения параметров:
a = 0.5 м; 0 = 2000 кг/м3; K0 = 14 ГПа; 0 = 8,4 ГПа; c1 = -0,09, c2 = 0,04 ГПа;
Su* = 36 МПа; 0 = 100 Па·с; = 1500 Па·с; C = 8·10-5(Па·с)-1; A = 1200 Па·с;
= 1,4; D* = 334,4 Дж/кг.
В расчетах было взято однородное распределение начальной поврежденности материала: 0 = 0(r,) = const. = 0,05. Для начального давления в порах p0 из (9) было взято значение (-)t=0, где (r,,t) - шаровая часть тензора напряжений.
Ниже представлены результаты расчетов с начальным распределением напряжений (12) для 1 = -75 МПа, 2 = -65МПа, 3 = -35 МПа; p = |3|.
На рис. 4 изображена зависимость скорости точек стенки скважины для разных значений угла от времени. Различие скорости стенки обусловлено тем, что 1 не равно 2.
Рис. 4. Радиальная скорость различных точек стенки скважины.
На рис. 5 изображены последовательные моменты времени для компонент тензора напряжений в зависимости от радиуса для различных углов при реализации первого подхода. Первый, второй и третий ряды – это радиальные, кольцевые и сдвиговые напряжения соответственно;
четвертый ряд – шаровая часть тензора напряжений. После снятия внутрискважинного давления стенка скважины устремляется к центу, и среда в радиальном направлении растягивается, и радиальные напряжения растут, оставаясь отрицательными (волна разгрузки бежит вглубь пласта). В угловом направлении среда наоборот сжимается, и кольцевые напряжения уменьшаются (в абсолютном значении увеличиваются). Сдвиговые напряжения для углов = 0 и =/2 равны нулю для всех моментов времени в силу симметрии задачи.
t = 0,00075 cек t = 0,00128 cек t = 0,00234 cек Рис. 5. Волновые картины для компонент тензора напряжений На рис. 6 изображено распределение удельной диссипации.
Критический уровень D* соответствует белому цвету. Видно, что область разрушения локализована в зоне, которая соответствует концентратору наибольшего сжимающего напряжения 2 (по вертикали).
зона разрушения Рис.6. Распределение удельной диссипации На рис. 7 представлена зависимость максимальной удельной диссипации в области в зависимости от времени.
Рис. 7. Максимальная удельная диссипация, Дж/кг Область разрушения при реализации второго подхода качественно меняет свой характер. Разрушение теперь осуществляется по типу магистральной криволинейной трещины (протяженные сравнительно узкие зоны разрушения), которые зарождаются на стенке скважины и прорастают вглубь пласта. На рис. 8 показаны распределения удельной диссипации энергии в последовательные моменты времени, по которым можно проследить развитие трещин – белый цвет соответствует критическому значению удельной диссипации. Количество трещин, направление и скорость их роста зависит от соотношения 1/2. Для симметричного случая (1=2) зона разрушения представляет собой кольцо вокруг скважины.
График зависимости глубины проникания трещины от времени для различных соотношений 1/2 представлен на рис. 9.
зарождения трещины поворота трещины Магистральная Рис. 9. Зависимость глубины проникания трещины от времени Зона разрушения при моделировании с учетом контактных взаимодействий представлена на рис. 10. Характер разрушения носит фрагментационный характер.
Зона разрушения Рис.10. Распределение удельной диссипации и фрагментационное 1. Впервые численно исследована задача необратимого динамического деформирования и разрушения горного пласта вблизи скважины в несимметричной двухмерной постановке с учетом как микроразрушения, так и с явным построением зон макроразрушения.
2. Сделаны выводы о характере и масштабах разрушений в пласте в зависимости от соотношения горных нагрузок вдали от скважины 1/2 и от степени реализации граничных условий на берегах макротрещин.
3. Показано, что упрощение алгоритма реализации граничных условий приводит к принципиально различным характерам разрушения пласта.
Поэтому, для получения физически реальных результатов, необходимо проведение численных исследований при максимально полной реализации граничных условий на берегах трещин, образующихся внутри пласта.
1. Kiselev A.B., Zacharov P.P. Computational simulation of irreversible deforming, micro- and macrofracture of rock in the vacinity of a borehole in its dynamical unloading // 8th World Congress on Computational Mechanics.
5th European Congress on Computational Methods for Coupled Problems in Science and Ingineering (Venice, Italy, 30 June – 4 July 2008). B.A.
Schrefler and U. Perego (Eds.).
Abstract
on CD-ROM. – Barcelona:
2. Захаров П.П., Киселев А.Б. К исследованию двумерной задачи Ломоносовские чтения. Научная конф. Секция механики. Апрель года. Тезисы докладов. - М.: Изд-во Моск. ун-та, 2008. - С. 98-99.
3. Захаров П.П., Киселев А.Б. Моделирование некоторых задач разгрузки в геомеханике // Современные проблемы газовой и волновой динамики.
Тезисы докладов международной конференции, посвященной памяти академика Х.А. Рахматулина в связи со 100-летием со дня рождения (21-23 апреля 2009, Москва, МГУ им. М.В. Ломоносова). - М.: Изд-во Моск. ун-та, 2009. – С. 39-40.
4. Захаров П.П., Киселев А.Б. Моделирование задач разгрузки в геомеханике // Ломоносовские чтения. Научная конф. Секция механики. Апрель 2009 года. Тезисы докладов. - М.: Изд-во Моск. унта, 2009. - С. 77.
5. Захаров П.П., Киселев А.Б. Численное моделирование динамики деформирования и разрушения горного пласта в призабойной зоне // Теория и практика расчета зданий, сооружений и элементов конструкций. Аналитические и численные методы. Сб. трудов международной научно-практической конференции / Московский государственный строительный университет. - М.: МГСУ, 2009. – С.
6. Kiselev A., Zacharov P. Irreversible Deforming, Micro- and Macrofracture of Rock in the Vicinity of a Borehole in its Dynamical Unloading // IV European Conf. On Computational Mechanics (Palas des Congress, Paris, France, May 16-21, 2010). – Abstracts on CD-ROM. – 1 p.
7. Киселев А.Б., Захаров П.П. Численное моделирование динамики деформирования и разрушения горного пласта в призабойной зоне // Ломоносовские чтения. Научная конф. Секция механики. Апрель года. Тезисы докладов. - М.: Изд-во Моск. ун-та, 2010. - С. 104.
8. Киселев А.Б., Захаров П.П. Численное моделирование динамики деформирования и разрушения горного пласта в призабойной зоне // Упругость и неупругость. Материалы Межд. научного симпозиума по проблемам механики деформируемых тел, посвященного 100-летию со дня рождения А.А. Ильюшина (Москва, 20-21 января 2011 г.) – М.:
Изд-во Моск. ун-та, 2011. - С. 476.
деформирования и разрушения горного пласта в призабойной зоне // Двойные технологии. – 2010. - №4. – С. 45-49.