Содержание к диссертации
Введение
1 Постановка задачи 12
1. Масштаб осреднении пороупругой среды. Модели Терцаги и Био 12
2. Хрупкое и квазихрупкое разрушение. Теория развития трещин Гриффитеа. Критерий разрушения Ирвина 17
3. Другие критерии разрушения. Критерии Кулоиа-Мора и Мизеса 27
4. Управляемое и неуправляемое развитие трещин. Обобщение теории хрупкого разрушения на случай пороупругой среды. 29
2 Численные методы 35
1. Метод опорных операторов 35
2. Решение систем линейных алгебраических уравнений 42
3. Многосеточный метод 45
3 Математическое моделирование задач пороупругости 57
1. Моделирование напряженно-деформированного состояния флюида при разработке месторождений нефти 57
2. Исследование задачи гидравлического разрыва 61
Заключение 69
Литература
- Хрупкое и квазихрупкое разрушение. Теория развития трещин Гриффитеа. Критерий разрушения Ирвина
- Другие критерии разрушения. Критерии Кулоиа-Мора и Мизеса
- Решение систем линейных алгебраических уравнений
- Исследование задачи гидравлического разрыва
Введение к работе
При разработке нефтяных и газовых месторождений используются различные технологии воздействия на нефтенесущие пласты, такие как закачка жидкости, газа и голевых смесей, а также щелочная обработка призабой-ной зоны, операции гидравлического разрыва пласта. Добыча углеводородов характеризуется интенсивным разбуриванисм и значительным движением жидкости в горной породе. Все вышеуказанное может приводить к возникновению зон повышенной деформации горной породы. Возникающие деформации приводят к срыву буровых колонн, движение жидкости приводит к нарушению равновесия и росту сейсмической активности. Для описания процессов движения флюидов и возникновения напряженно - деформированного состояния используются различные математические модели. Эти модели очевидным образом должны включать уравнения фильтрации и уравнения теории упругости, при этом тензор напряжений должен содержать часть, зависящую от порового давления. Предложенное К, Терцаги соотношение, содержащее напряжение скелета, не является содержательным, так как смещения скелета не наблюдаемы. Уравнения движения флюида в пористой среде может быть описано как движение многофазного континуума. При этом записываются уравнения непрерывности, уравнения сохранения импульса и энергии. Силы, действующие на твердую фазу - это упругие силы и силы межфазного взаимодействия [Рахматулин Х.А., 195G). На флюид действует сила давления и также силы межфазного взаимодействия. Полные напряжения определяются согласно формуле Терцаги [Tcrzaghi К., 1968]. Вышеизложенная модель позволяет описать эволюцию многофазного континуума при больших скоростях относительного движения, однако она содержит ряд величин, которые должны быть доопределены, как то, сила межфазного трения и потоки тепла по фазам. Если скорости движения как скелета, так и флюида малы, то можно построить модель, в которой скорость флюида выражается законом Дарси, а тензор напряжений определяется не через тензор деформаций ске-
лета, а через тензор деформаций в целом [Biot М.А., 1962]. При быстром нагружении трещиноватых пороупругих слоев в них возникает дилатанси-онный э(])фект, он сопровождается резким понижением норового давления. При этом могут образовываться крупные трещины и разломы. Уравнения описывающие данные эффекты были получены в работе [Райе Дж., 1982] и в дальнейшем развиты в работах [Rice J.R., Cleary J.M., 1976] и др..
В настоящее время большинство нефтегазовых месторождений России находятся па стадии завершающей разработки. Такие месторождения характеризуются общим истощением и большой обводненностью нефтедобывающих скважин. Новые месторождения, вовлекаемые в разработку, содержат трудно извлекаемые запасы углеводородов. Коллекторы характеризуются низкой проницаемостью и слабым дренированием. В течение длительного срока эксплуатации скважин, параметры призабойпой зоны значительно ухудшаются. Это связанно изменением проницаемости, выпадением парафинов и асфальтепов и значительной обводненностью скважин. Эффективность работы таких скважин за время эксплуатации значительно уменьшается. Одним из основных методов интенсификации разработки сложных и проблемных нефтегазовых месторождений является гидравлический разрыв пласта (ГРП). Метод основывается на механическом воздействии на горную породу. Под действием избыточного давления закачиваемой жидкости порода разрывается по поверхности минимальной прочности. Для этого в породу закачивается жидкость с расходом, много превышающим способность породы к поглощению жидкости. После разрыва пласта образуется трещина или система трещин, обладающие высокой проницаемостью. Зона растрескивания повышает зону дренирования в зоне работы скважин, тем самым достигается высокая эффективность добычи. В настоящее время около трети запасов углеводородов можно извлечь только с использованием этой технологии [Константинов С.В.,Гусев В.И., 1985].
При изучении процессов гидравлического разрыва пласта выделяют два направления:
Эвристические модели
Геомеханические модели
Для проведения расчетов в настоящее время преимущественно используются эвристические методы. В первую очередь, это связано со слабым развитием геомеханических моделей и сложностью решения соответствующих
задач. Успех такого подхода базируется на большом количестве экспериментального материала и простоте использования метода. Однако, эвристические модели не имеют строгого научного обоснования, основываются на статистике, что является минусом данного подхода. Попытка построить геомеханическую модель была предпринята в работе |Желтов Ю.П., 1955]. Задача рассматривалась в плоской постановке в упругом приближении в изотропной и однородной среде. Изучалась трещина, в которую из скважины нагнетается вязкая жидкость (смесь геля и пропапта). Движение жидкости описывается уравнениями Стокса в приближении смазки. В такой постановке давление меняется только вдоль трещины, поперек трещины давление считается постоянным. Основные предположения сделанные в работе [Жслтов Ю.П., Христиапоиич С.Л., 1955] -это принцип подобия трещины и условие плавного смыкания берегов трещины. Данное обстоятельство значительно облегчает решение задачи. В такой постановке, при известном постоянном давлении на берегах трещины в однородной и изотропной среде, может быть получено асимптотическое решение. Недостатком данного подхода является постулирование формы трещины - при увеличении размеров трещина не меняет своей формы. Впоследствии модель была усовершенствована в работах Г.И. Барен-блатта [Баренблатт Г.И., 19бі|, в первую очередь, введением сил сцепления вблизи концов трещины. В дальнейшем модель подвергалась изменениям, в настоящее время задача гидроразрыва выглядит следующим образом [Kovscek A.R., Johnston R.M., 199G|. Рассматривается пороупругая изотропная однородная среда с единичной трещиной. В трещине движется вязкая жидкость, описываемая уравнениями Стокса или Навье-Стокса, жидкость внутри трещины гидроразрыва предполагается не ньютоновской [Valko P., Economides М.Л., 1995].
Давление на границе трещины предполагается равным нормальным напряжениям вдоль трещины (с обратным знаком). Задача разбивается на два этапа. На первом этапе решается прямая задача нахождения зависимости давления от потока в скважине (уравнения Стокса или Навье-Стокса). На втором этапе решается обратная задача нахождения размеров трещины при известных закономерностях зависимости давления от потока. Недостатки данного подхода очевидны. Необходимо решать две задачи - прямую и обратную. Возникают определенные сложности с решением обратных задач в связи с их неустойчивостью. Кроме того, существенным недостатком данного подхода является то, что модель не учитывает филь-
трациоппые движения жидкости вне трещины гидроразрыва. Кроме того данная модель не является самосогласованной, а именно, на границе трещины и окружающей среды предполагается выполнение условия равенства нормального напряжения и давления жидкости, но не предполагается равенства смещений величине раскрытия трещины.
Для устранения данных недостатков возможно использование более пол-пых уравнений, описывающих пороупругую среду (уравнения Био). Использование пороупругих уравнений позволяет правильно и согласованно описывать движение породы, учитывая его фильтрационные и упругие характеристики, что затруднительно в рамках других моделей. При использовании такого подхода, зная распределение давления внутри трещины (посредством решения уравнений Стокса), решается полная связанная задача Био. В этой постановке не постулируется ни форма трещины, ни ее размеры. Размеры и форма трещины будут зависеть от возникающих полей напряжений и фильтрационных потоков. Для определения раскрытия трещины и зон разрушения породы используются различные критерии разрушения. Практика показывает, что после проведения операции ГРП не наблюдается одна большая трещина гидроразрыва, образуются значительные зоны направленной трещиноватости [Курьяпов Ю.А. и др., 2001]. Таким образом, в результате проведения гидроразрыва образуется аномальная, неоднородная зона повышенной трещиноватости. Поэтому, используя критерии разрушения, правильнее говорить не о трещине гидроразрыва, а о зоне искусственной трещиноватости. Соответственно, используя критерии разруигения, необходимо определять такие зоны. К преимуществам модели можно отнести то, что в данной постановке возможно учесть анизотропию породы и ее неоднородность. Недостатком подхода является большие вычислительные затраты. Однако, в свете использования мощных вычислительных ресурсов, таких как кластерные системы, эта проблема преодолима. В некоторых специальных случаях можно пренебречь влиянием напряженного состояния на движение флюида, в этом случае размерность задачи меньше исходной, и трудоемкость проведения расчетов снижается. В такой постановке уравнения пороунругости поддаются расщеплению на фильтрационную и упругую части, тогда исходные уравнения значительно упрощаются [Каракин А.В. др., 2003]. В настоящей работе рассматривается задача пороунругости (полная связанная задача Био) в результате чего определяется распределение порового давления и полей напряжений. С помощью критериев разрушения определяется зона возможного разрушения,
что позволяет описать возможную область искусственной трещиноватости. Среда характеризуется значительной неоднородностью физических параметров, таких как проницаемость и пористость. Земная кора имеет сложную геометрическую структуру, пронизана трещинами и разломами. Для моделирования таких сложных сред необходимо использование адаптированных к среде сеток или специальных высокоточных алгоритмов усреднения. В связи с этим, к алгоритмам предъявляется требование корректного описания разрывов. Существенной особенностью проблемы построения алгоритмов в средах сложной структуры является обеспечение непрерывности потоков рассматриваемых неличин. В многочисленных работах отечественных и зарубежных авторов эта задача решается посредством метода конечных объемов, при этом сетка предполагается, как правило, адаптированной к структуре среды, а значения величин относят к ячейкам. При определении значений потоков величины коэффициентов вычисляются как среднее арифметическое для компонентов потока параллельных границе или как среднее гармоническое, если направление, вдоль которого определяется значение потока, пересекает границ}' разрыва коэффициентов. Трудности построения таких алгоритмов особенно ярко проявляются на не ортогональных сетках. При этом имеет место значительная ориентаци-онная ошибка. Для преодоления данной проблемы в методе конечных объемов используются различные весовые функции, компенсирующие ошибки, возникающие при прямой аппроксимации потоков, т.е. для аппроксимации используются адаптированные модифицированные алгоритмы, основанные на методе конечных объемов [Колдоба А.В. и др, 1985]. В средах со сложной геометрией трудно обойтись регулярными сетками. Геологическая структура пород такова, что нередки выклинивания пластов. В таких случаях использование нерегулярных сеток становится необходимостью. Использование алгоритмов, основанных на методе конечных объемов, в этом случае затруднительно или невозможно. Техника построения разностных схем, преодолевающих данные препятствия, была разработана А.А. Самарским и его сотрудниками [Самарский А.А. и др., 1996]. Основная идея метода опорных операторов - это согласованное в силу теоремы Гаусса-Остроградского определение разностных аналогов операторов градиента и дивергенции. Метод опорных операторов на разрывах имеет аппроксимацию порядка 0(h), в то время как обычные методы конечных объемов только 0(1). До последнего времени в работах, посвященных опорным операторам, рассматривались только среды с непрерывными коэффициентами
в ввиду сложности аналитического исследования алгоритма. Б настоящее время проведено качественное изучение алгоритма в случае анизотропных и неоднородных, в том числе с разрывными коэффициентами, сред [В.П. Мясников, М.Ю. Заславский, А.Х. Пергамент, 2004].
Для обеспечения точности расчетов часто используется большое количество узлов расчетной сетки. При подготовке инженерных расчетов становится актуальной задача минимизации времени, затраченного на расчет, поэтому численные алгоритмы и методы должны обладать достаточной скоростью сходимости для выполнения данных требований. При построении неявных схем в задачах математической физики возникает необходимость решения большой системы линейных алгебраических уравнений. Данные задачи характеризуются жесткой матрицей и плохой обусловленностью. В связи с этим, для решения таких систем должны быть использованы специальные методы. Для решения задач с симметричными матрицами использованы хорошо зарекомендовавшие методы, такие как метод сопряженных градиентов или верхней релаксации. Общеизвестно, что метод конечных объемов в случае пеортогоиальных сеток не обеспечивает симметричность полученной матрицы. В случае несимметричной матрицы выбор метода затруднителен. Для такого класса матриц большинство итерационных методов не эффективны и имеют крайне малую скорость сходимости. В связи с этим, создание и использование эффективных методов является актуальной задачей. Одним из эффективных методов уменьшения времени обращения матриц является использование предобуславли-вателей. Основная идея метода заключается is том, что исходная система линейный алгебраических уравнений трансформируется в другую систему, характеристики которой значительно лучше исходной, тем самым итерационный метод сходится значительно быстрее [Saad Y., 1995]. В настоящее время данное направление бурно развивается, в свободном доступе существуют библиотеки, позволяющие использовать различные предобуслав-ливатели и итерационные методы (Sparskit, Aztec, Sparslib). В работе для решения системы линейных алгебраических уравнений использованы пре-добуславливатели ILUO, ILUK, ILUTP, в качестве итерационного метода использовался метод минимальных невязок.
Принципиально иной подход для уменьшения времени счета и трудоемкости эллиптических задач был использован Федоренко Р.П.. В основе разработанного многосеточпого метода [Федоренко Р.П., 1994] лежит идея использования многоуровневых сеток. Известно, что несколько простых
итераций или итераций по Зсйдслю эффективно погашают высокочастотную составляющую ошибки и сглаживает невязку. Поэтому формируется задача вычисления поправки па более грубой сетке, так как поправка стала более гладкой функцией. Полученные результаты на различных уровнях сеток интерполируются на исходную расчетную сетку и из них составляется решение задачи. Скорость сходимости многосеточного метода определяется выбором уровней сеток, количеством итераций и наличием хорошей интерполирующей функции. Несмотря на большую трудоемкость метода, особенно в случае не ортогональных сеток, метод является одним из самых передовых и часто используемых.
Цель и задачи. Целью данной работы является изучение связанной и не связанной задач Био, описывающих пороупругую среду, и выявление особенностей траектории гидравлического разрыва в неоднородной среде. Построение моделей взаимодействия основной трещины гидравлического разрыва пласта с геологическими нсодпородностями. Разработка и реализация эффективных численных методов на основе метода опорных операторов и мпогоесгочиого метода для решения уравнений Био.
Методы исследований. Построен эффективный алгоритм решения системы уравнений пороупругости. Для изучения процессов реализованы трехмерные программы решения связанной и не связанной задачи Био. Для случая не связанной задачи применяется многоссточный метод решения уравнений фильтрации.
Для исследования процессов гидравлического разрыва пласта реализована двумерная программа. Для определения возможных зон разрушения используются критерии разрушения.
Научная новизна. Научная новизна определяется решением трехмерной связанной задачи Био. Рассмотрены методы построения алгоритмов в среде с разрывными коэффициентами. Исследована задача взаимодействия трещины гидравлического разрыва пласта с естественными геологическими нсодпородностями. Показано изменение траектории трещины гидроразрыва в неоднородной среде.
Апробация работы. Оснотшые результаты работы докладывались па международной конференции „XXII Европейский симпозиум по повышению нефтеотдачи пластов", Казань 2003; на XV Всероссийской конференции "Теоретические основы и конструирование численных алгоритмов для решения задач математической физики с приложением к многопроцессорным системам ", посвященной памяти К.И. Бабеико; на научных семинарах
им. Бабенко К.И. ИПМ РАН.
Публикации. По теме диссертации опубликовано пять печатных работ, в том числе две в соавторстве. Из них 1 статья в российском журнале, 2 статьи в препринтах ИПМ им. М.В. Келдыша РАН, 2 работы в сборниках трудов докладов международных и всероссийских конференций.
СОДЕРЖАНИЕ РАБОТЫ
Первая глава работы посвящена математическому описанию поро-упругой среды. Приведено построение математической модели па осново уравнений Био. Уравнения, описывающие пороупругую среду, включают в себя систему уравнений упругости и фильтрации, при этом фильтрационные качества породы зависят от шаровой части тензора деформаций. В свою очередь напряжения и деформация определяются градиентом поро-вого давления.
Рассмотрены задачи теории трещин. В однородной, изотропной средо для одиночной трещины могут быть получены асимптотические решения. В случае пороупругой однородной и изотропной среды уравнения системы Био расщепляются па несвязанные уравнения теории упругости и фильтрации, тем самым для задачи Био, в такой постановке также, могут быть получены асимптотические решения.
Во второй главе диссертации рассмотрены численные методы, используемые для решения задач пороупругости и фильтрации. Алгоритмы основаны па методе опорных операторов и многосоточпом методе. Реализованные программы адаптированы на случай неоднородной среды и используют неортогональные расчетные сетки.
Проведено исследование сходимости многосеточного метода. Скорость сходимости метода определяется в сравнении с методом верхней релаксации, имеющим тот же порядок сходимости. При интерполяции данных с различных уровней сеток используется полилинейная интерполяция. В случае ортогональной сетки может быть использована интерполяция Кинга.
Третья глава диссертации посвящена вопросам изучения пороупругих эффектов. Реализованная программа оттестированы на известных, классических задачах теории упругости и теории трещин. Проведено моделирование трехмерного коллектора с разломом. Задача решается в геологическом кубе, куб содержит в себе нефтенссущий коллектор с разломом. Фильтра-
ционпые потоки и значения давления определяются в зоне фильтрации, упругие поля напряжений и деформаций во всем кубе. Показаны графики распределения полей напряжений.
Рассмотрена двумерная задача взаимодействия трещины искусственного гидравлического разрыва с геологической неоднородностью (трещиной). Взаимодействие осуществляется иод различными углами и при различных расстояниях между трещинами. Приведены соответствующие графики. Получены пространственные распределения компонент тензора напряжений и функция давления. Для определения возможных зон разрушения использован метод Кулона-Мора. Приведены графики, демонстрирующие различные сценарии разрушений в зависимости от типа взаимодействия трещин.
Хрупкое и квазихрупкое разрушение. Теория развития трещин Гриффитеа. Критерий разрушения Ирвина
Большой вклад в изучение процессов трещинообразования внесли Гриффите, Ирвин, Орован. Их труды в первую очередь относятся к хрупкому и квазихрупкому разрушению. Под хрупким разрушением обычно понимается разрушение при котором после разрушения части тела можно сложить, так чтобы составленное тело совпадало с исходным [Седов Л.И., 1950]. Для природных материалов более характерно квазихрупкое разрушение, при котором образуется малая пластическая зона. Эксперименты проведенные Ирвином и Орованом показали, что в большинстве практически важных случаев разрушение происходит квазихрупким образом, что позволяет использовать теорию малых деформаций и хрупкого разрушения для корректного описания процессов трещипообразоваиия. Далее рассматриваются трещины и щели уже имеющиеся в материале, вопросы начального возникновения трещин не рассматриваются.
Определяющее значение в изучении процессов образования трещин имеют вопросы касающиеся критериев разрушения. Первые работы по теории трещин принадлежат Гриффитсу. Гриффите [Griffith А,А., 1920] предложил энергетический критерий разрушения трещин, при этом Гриффите исходил из решения об эллиптическом вырезе в упругой неограниченной плоскости. Решение получено устремлением полуоси эллипса к нулю.
Данное решение поддается обобщению на случай эллиптического и прямолинейного выреза выбором соответствующего отображения [Седов Л.И., 1950]. Действительно, рассмотрим функцию Жуковского: 0 т [1.27)
Окружности ; = 1 соответствует эллипс с центром в начале координат и полуосями а = Я(1 + т), Ь = Н{1 — т). Подобрав соответствующие значения a, b можно получить любой эллипс, в предельном случае т = 1 эллипс обращается в двусторонний отрезок, заключенный между точками ±2R. Применив соответствующее обратное отображение можно свести задачу с эллиптическим вырезом к задаче о круге, подробно данные вопросы изучены, например, в работах [Седов Л.И., 1950, [Мусхелишвилли Н.И., I960]. В случае всестороннего растяжения бесконечно тонкого эллиптического выреза длиной 2а, распределение нормальных напряжений
Здесь Е, v коэффициенты Юнга и Пуассона соответственно. Согласно решению, при приближении к носику трещины, нормальные напряжения бесконечно возрастают (рис.1.2), при этом берега трещины совершают вертикальное перемещение определяемое формулой: блем прочности и разрушения Гриффите основывался па универсальных уравнения термодинамики. Закон сохранения энергии для тела в конечных размеров в обідом случае имеет вид: dE + dU = dA + dQ-\ dQ (1.30)
Здесь Е- кинетическая энергия, U - полная внутренняя энергия, dA -работа объемных и макроскопических сил, dQ - внешний приток тепла, dQ - „особый" приток энергии, обусловленный, например, химическими реакциями. Для развивающейся трещины в закон сохранения (1.30) необходимо добавить член dA - энергия в особых точках, возникающая за счет движения краев трещины, закон сохранения примет вид: dE + dU = dA + dQ + dQ + dA" {1.31) Далее, в предположении хрупкого и квазихрупкого разрушения, считается, что величины dE, dQ, dA не изменяются, тогда закон сохранения примет вид: dU = dA" + dQ (1.32)
Для расширения трещины необходимо преодолеть силы сцепления, т.е. совершить некоторую работу. Для образования единицы свобоной поверхности необходимо совершить работу 7, величина 7 называют плотностью поверхностной энергии, ее можно считать константой в данных условиях. Соотвотствсіпю при раскрытии трещины длиной 1а выделяется энергия А — А уа Согласно [Griffith А.А., 1920], напряжения совершают работу W: W = 2 (-±P0juA =-(1- J) 3L (х.зз) —а Если трещина растет, то освобождаемая энергия А расходуется на образование новой поверхности: dW = dA (1.34) соотношение можно переписать в виде: -?-(W + A) = 0 (1.35) Учитывая количество энергии необходимое для раскрытия трещины и работу совершаемую напряжениями окончательно получим формулу: 7га(1 — Vі)
Эта формула при заданной длине трещины а определяет критическое значение напряжений Р на границе трещины, при котором будет происходить рост трещины. Формулу впервые получил Гриффите в 1920 году. Рассмотрим упругое тело П, ослабленное трещиной О, с поверхностью Е, на границе области Q, приложены нагрузки Р{0), для каждой точки пространства определены смещения щ и деформации Cfj. Под действием нагрузок внутренние точки совершают перемещения Агц, при этом поверхность трещины увеличивается на величину dH. В этом случае, измепие 0 Я1 ґ0 6L а) Ь) Рис. 1.3: Изменение поверхности трещины внутренней энергии определяется формулой [Качаної! Л.М., 1974]: AU= I P{S)AuidSi (1.37)
В случае симметричной трещины поверхность Е — 2+ + Е_, Е+ — _, где +, _ -верхний и нижний берег трещины соответственно. Если приложенные силы Р симметричны относительно оси трещины : AU = IJ P{S)AuidSi ША (1.38) В 1957 году в работе [Irwin G.R., 1957] Ирвин интеграл в правой части представил в виде: AU = і ( P(S)AUidSi - -QdZ+ (1.39) JE_, в предположении существования некоторой постоянной величины Q, постоянной для данной трещины в заданных условиях. Множитель Q имеет структуру [Качанов Л.М., 1974], [Броек Д., 1973],: Q = ск2 (1.40) где к- коэффициент интенсивности напряжений. В силу того, что изменение внутренней энергии может быть определено по формуле dU — 2rydH,i окончательно приходим к формуле Ирвина, для определения критического значения коэффициента интенсивности напряжений к#: k = J- = k (1.41)
Трещина начинает распространятся если коэффициент интенсивности напряжений достигает критического значения & . Необходимо заметить, в этом случае критерий разрушения носит локальный характер, поскольку оно зависит от полей напряжений вблизи вершины трещины. Таким образом энергетический (глобальный) критерий Гриффитса можно заменить силовым (локальным) критерием вблизи вершины трещины. Позднее Ирвин показал эквивалентность обоих критериев. Необходимо заметить, что практические исследования не подтверждают того что к строгая константа материала, в общем случае коэффициент может зависеть от многих условий. Однако, критерий Ирвина существенно лучше других параметров материала, характеризующих чувствительность разрушений к концентрации напряжений.
Другие критерии разрушения. Критерии Кулоиа-Мора и Мизеса
В настоящее время требования к численным решениям существенно ужесточились. В первую очередь это связано с широким внедрением вычислительных методов для сложных технологических расчетов. Увеличение точности за счет увеличения размерности сетки не всегда оправдьшет себя. Во-первых увеличение размерности сетки не гарантирует корректное описание геометрии. Кроме того, в процессе расчета область существенно может менять свою геометрию, например при наличии сильных сдвиговых деформаций. В этом случае не редко сетки разваливаются с потерей аппроксимации, что приводит к срыву расчета. В связи е возникающими такого рода проблемами, в практике широко используются псоргогональпые и нерегулярные расчетные сетки адаптированные с структуре среды. Аппроксимация уравнений на таких сетках достаточно трудная задача, раз-постная схема, получаемая в ходе аппроксимации, должна обладать рядом необходимых свойств для применения численных методов решения. Такими свойствами являются консервативность схемы, ее монотонность и положительная определенность.
Техника построения разностных схем для уравнений математической физики, обладающие необходимыми свойствами, были развиты в работах группы А. А. Самарского. В качестве базового в работе выбран метод опорных операторов. Основная идея метода опорных операторов это согласованное определение разностных аналогов операторов дивергенции и градиента.
Рассмотрим метод опорных операторов на примере скалярного эллиптического уравнения. div(-fcgrad(u)) = / (2.1) Для согласованного определения операторов, один оператор аппроксимируется непосредственно, второй таким образом чтобы удовлетворить разностному аналогу основного интегрального тождества: vdivqdV= (v{pdS)+ f(q,p)dV (2.2) О 80 О где q = — kgrad(u) Пусть функции U, V-сеточные функции определенные
Расчетная сетка в узлах расчетной четырехугольной сетки М. Сетка состоит из множества ячеек П, узлов и соединенных множеством ребер 7- Предполагается постоянство коэффициента к внутри ячейки. Общий вариант, позволяющий учитывать наличие разрывов коэффициентов, разработан в работе [В.П. Мясников, М.Ю. Заславский, А.Х. Пергамент, 2004]. На ребре а Є -у определим величины Qa, Ра: Qa = -ДС/ - -U(B) + U{A) Pa = -AU = V{B) + V{A) (2.3) То есть Qa, Ра есть разностный аналог величин qa = Jgradu и ра = а J gradu, соответственно. Таким образом, задача свелась к аппроксимации а второго интеграла в уравнении (2.2). Этот интеграл представляет собой билинейную форму величин Qa, Ра и, следовательно, может быть представлен в виде: J2gabQaPb = (Q,P) (2.4) 0ЄЇЗ где раЬ-матрица билинейной формы для ячейки G, здесь сумма по всем ребрам, ограничивающим ячейку G. В случае если на границе области V\QO—О интеграл по границе J v(pdS) также равен нулю тогда разностная аппрок симация интегрального тождества (2.2) примет вид: [ vdivqdV = J2Vi (div (kSrad(u)))i vi = Yl i7 u (2-5) JQ гЄи і
В правой части уравнения слаї-аемое aaQa суммируется по всем ребрам во всех ячейках содержащих ребро а. Здесь аа 1 если а выходит из узла г, соответственно аа = — 1 если а входит. Величина Qa определяется соотношением [А.А. Самарский и др., 1996]: Qa = GQ = gabQb, (2.6) здесь суммирование производится по всем ребрам b имеющих общую вершину с ребром а. Оператор G порожден семейством матриц даЪ} определенных для каждой ячейки в, и действует на всей области решения О, Таким образом разностная задача для скалярного эллиптического уравнения свелась к виду: oaQa = fi= I fdV (2.7) Здесь Vi некоторой присоединенный объем к узлу і, на рисунке выделен пунктиром. Определение присоединенного объема обсуждается далее.
Точность любого метода определяется точностью аппроксимации интеграла энергии и близостью потоков в соответствующей норме. Для достижения желаемой аппроксимации необходимо правильно построить метрический оператор даЬ и выбрать присоединенные объемы V%. В случае постоянных коэффициентов к во всей области способ построения метрического оператора предложен в работе [А.А. Самарский и др., 1996], построение даЬ. В случай разрывных коэффициентов метод построения метрического оператора предложен в работе [В.П. Мясников и др., 2004]. Интеграл J kgrad(u)gvad(v)dV, являющийся энергетическим скалярным произве-о дением и и v, можно представить в виде (2.4). Аппроксимация интеграла энергии в шотнетст вую ідей норме обеспечивает сходимость метода в слабом смысле. Для достижения сильной сходимости необходимо выполнения условия равенства потоков. Для метода опорных операторов и двумерном случае удается доказать сильную сходимость. В трехмерном случае этого достичь не удается и алгоритм имеют только слабую сходимость. Необходимо заметить, что в трехмерном случае обеспечивается консервативность схемы и функция имеет аппроксимацию порядка 0(h), потоки аппроксимируются величиной порядка 0(1). Для определения вида метрического оператора даЬ потребуем точного совпадения непрерывной энергии /fc(gradu) У и разностной д QaQb на функциях линейной оболочки е L(0) = зрап(1,х,у). В этом случае [В.П. Мясников. М.Ю. Заславский, А.Х. Пергамент, 2004] действительна теорема. Теорема 1. Существуют такие дпЬ такие, что J к (graclu) dV = gabgaQb є для и Є L(&), и они могут быть представлены в виде: даЬ = к Зф(і а,і ь) (2-8) Ф
Суммирование производится по всем узлам ячейки G. Выражение (la,l j-евклидово скалярное произведение векторов, сопряженных к векторам 1а, 1ь с общей вершиной ф. Вектор 1а по модулю равен длине ребра а и направлен вдоль положительного направления (рис.2.2). S$ присоединенные объемы к вершинам ф ячейки 0: = 5@, / -значение коэффициента в ячейке
Как уже говорилось, тождественное выполнение равенств непрерывной и разностной энергии влечет слабую сходимость, специальным подбором присоединенных объемов Бф и Vi можно обеспечить сильную сходимость. Здесь необходимо пояснить разницу между объемами 5 и Vi. Проиллюстрируем это графически (рис.2.3). Область выделенная серым цветом присоединенный объем к узлу А, состоит из Бф принадлежащих различным ячейкам, тем самым Уд есть некоторой объем новой ячейки, для которой фактически пишется баланс.
Решение систем линейных алгебраических уравнений
В последние годы широкое применение при решении задач математической физики получил мпогосеточный метод (MULTIGRID). Метод является асимптотически самым быстрым по сравнению со всеми приведенными выше методами. Так, для уменьшения нормы погрешности в е раз требуется всего СМ арифметических операций. Константа С не зависит от количества неизвестных М системы Аи — Ь.
Идея метода заключается в использовании нескольких уровней сеток. На исходной сетке ищется решение, далее производится коррекция решения на вспомогательной, более крупной, сетке. Поясним метод па примере одномерного уравнения Пуассона. ( Д" = ё = /(г) (2.2С) I U\ y — ф 4G Запишем разностное уравнение: ик-1 - 2щ + ик+1 = Ik (2.27) /г2 Метод установления примет вид: 4+1 = 4 + г(Дл4 - Л) (2-28) Далее введем фундаментальные объекты - погрешность Vk и невязку гЦ.: 4 = «І-«ь гІ = ДА«І-Л (2.29) t4 Httt- приближенное и точное решение задачи (2.26) вузлах расчетной сетки. Для нахождения погрешности vj: получим уравнения: { = (2.30)
Как отмечалось выше, итерационный процесс (2.28) сходится очень медленно, однако, он эффективно погашает высокочастотные компоненты и сглаживает невязку [Федоренко Р.П., 1994]. Тем самым функция v после нескольких итераций становится гладкой сеточной функцией.
Введем операторную запись уравнения. Определим, оператор D, переводящий сеточную функцию в такую же функцию (Аи)г ,тогда для сеточных функций получим операторную запись v1+1 = (Е + TD)V\ Будем использовать следующие факты из теории дискретных рядов Фурье. Всякая функция Vk, рапная на границе нулю, имеет представление: где z собственные функции разностного оператора А, при этом = Х )2} (2-32) і „ і 1/2 Ml Теперь разлагая в ряд Фурье начальную погрешность v получим: v1 = (E + TD)V = {Е + TD) CpzP = 22 СР(Е + TD)ZP = = ?(!\ )Z? P (2 33) v Далее воспользуемся тем, что собственными функциями одномерного опе ратора D являются функции s = sin крп . Им соответствуют собствен ные значения Хр = T sin (—гт ] [Федорснко Р.П., 1994 . (Iі \2N / Тогда для погрешности получим формулу; . 2 / ртг ІЇ \Ш V = Е -4 sin кртг Построим график эволюции погрешности. Для определенности N = 10. N Условно разделим собственные функции па две части: гладкие р — N и негладкие р —, т.е. условно разделим низкие и высокие гармоники [Федоренко Р.П., 1994]. На графике (рис.2.5) приведены два графика за 0.6- —/ "X. 0.5- і і і! 0.4- 1 \ 0.3 і ііІ \\ 0.2- 1 іі іі \ 0.1-0 V і Т і і т2 4 р 6 8 Рис. 2.5: Эволюция погрешности виснмости погрешности от номера гармоники для различного количества итераций, на графике справа количество итераций меньше. Из графиков видно, что с увеличением числа итераций вклад высокочастотной составляющей погрешности эффективно погашается итерационным процессом. На остальной части спектра компонента погрешности убывает достаточно медленно. В целом такой итерационный процесс неэффективен. Однако небольшое число таких итераций "сглаживает невязку". Поэтому для определения v можно использонать более грубую расчетную сетку. Новая за 48
дача (2.30) для определения погрешности v содержит меньшее количество неизвестных и число обусловленности для системы меньше, соответственно она проще задачи (2,28).
По аналогии с вышеизложенным, в задаче (2.30) можно определить погрешность и невязку и получить новую задачу, которую можно решать на еще более крупной сетке. Таким образом, получим цепочку задач на нескольких сетках. Интерполируя результаты на основную сетку, получим решение исходной задачи (2.2G).
Изложенная схема относится к задаче с постоянными коэффициентами или, в крайнем случае, с непрерывными коэффициентами. Если коэффициенты разрывные, то по вышеизложенной схеме можно объединять только ячейки в области непрерывности коэффициентов. Действительно, при реализации метода неоднократно производится интерполирование, например, с использованием полилинейных функций. Подобное интерполирование является корректным только в том случае, когда функция достаточно гладкая. Поэтому в нижеизложенной схеме метода Multigrid для разрывных коэффициентов ячейки объединяются таким образом, чтобы не содержать внутри линий (или поверхностей) разрывов.
Алгоритм метода допускает обобщение, при усложнении задачи за счет неоднородности коэффициентов, произвольности вида области, применения неструктурированных сеток [Фсдоронко Р.П., 1994]. Для исследования многосеточного метода рассмотрена задача: div(Q) =f(x,y,z) Q=-k(x,y,z)gmdU t2 34) Область решения 3 может иметь сложную геометрическую структуру, коэффициент к- кусочно-непрерывная функция. В такой ситуации целесообразно использовать сетки, адаптированные к структуре разрывов коэффициентов, причем значения функций относятся к ячейкам, а потоки вычисляются на гранях. Для аппроксимации уравнений в этом случае используется модификация метода конечных объемов. Рассмотрим аппроксимацию подробнее. Введем семейство расчетных сеток следующим образом. Пусть в области решения построена сетка G, состоящая их ячеек Q, к которым отнесены сеточные функции. Множество вершин, т.е. узлов сетки, обозначим через ш. Пусть па гранях а определены значения нормальной компоненты потоков q.
Исследование задачи гидравлического разрыва
При изучении процессов гидравлического разрыва пласта, в некоторых случаях, задачу рассматривают в упругом приближении. Считается , что основная масса жидкости, закачанной в трещину гидроразрыва, расходуется на работу по увеличению трещины, при этом фильтрационными процессами пренебрегают. В более сложной постановке учитывается потеря жидкости гидроразрыва за счет фильтрации, но процессы фильтрации при этом, также не рассматриваются. Задавая геометрические размеры трещины и потери жидкости гидроразрыва за счет фильтрационных процессов можно определить усилия на берегах трещины. В этом случае возникающие поля напряжений и смещений хорошо изучены.
Реализованный алгоритм протестирован па классических решениях теории трещин в упругой постановке. Рассмотрим бесконечную упругую пластину ослабленную вырезом. Пусть к берегам трещины приложены нагрузки Р, известно, что в зависимости от вида нагрузок развитие трещин может проходить по двум сценариям. В первом случае для роста трещины необходимо увеличивать нагрузки приложенные к берегам трещины. Во втором случае рост трещины лавинообразен, при этом при увеличение длины трещины нагрузки на берегах трещины не возрастают. На графиках (рис.3.4), (рис.3.5) продемонстрировано сравнение численного и асимптотического решения. В случае одиночной трещины, для роста трещины необходимо увеличивать усилия на берегах Р, в то время как для случая неустойчивого развития трещины, необходимые усилия убывают с ростом длины трещины а. Соответственно, при устойчивом развитии трещины Р возрастает с ростом а. На графиках видно хорошее совпадение асимптотического и численного решения.
Важной задачей при проведении гидравлического разрыва пласта, является определение зон разрушения и направление эволюции трещины. Задача особенно актуальна в средах сложной геометрии, содержащих массу неоднородностей и особенностей, каковыми без сомнения являются природные коллектора. В таких случаях на процесс формирования трещины гидроразрыва большое влияние оказывают природные разломы и трещины, поэтому при моделировании необходимо учитывать возможность такого влияния. В этом случае возможности предыдущего метода сильно ограничены, в силу большого влияния процессов фильтрации на формирование трещины. Для моделирования процессов связанных с гидравлическим раз Г — 0.б, проницаемость к = U.UU5 Дарси, а для разлома: А = 0.2 10 Па, ц = 0.7- 109Па, ; = 0.95, к = 0.1 Дарси. В трещину закачивается жидкость разрыва. На границе трещины в качество граничных условий используем характерное распределение давлений Р(х,у). Напряжения будем отчитывать от состояния гсостатического равновесия. Ниже представлена схема взаимодействия трещин (рис.3.6).
При моделировании используются многоблочные сетки, с целью построения адаптивных сеток. Такая сетка строится в зоне разлома, во всей остальной области сетка не адаптивна. У трещины гидроразрыва адаптивность сетки но требуется, так как среда предполагается однородной. Вблизи зон концентрации напряжений сетка имеет сгущающийся характер (рис.3.7). При моделировании использовались различные размерности сеток, максимальная размерность составила 3 104 узлов.
При решении задачи о напряженно - деформированном состоянии одиночной трещины используется следующая схема. Задача условно разделена на два этапа. На первом этапе решается система уравнений теории поро-у пру гости, на каждом временном шаге определяются и значения смещений и поле напряжений. Для этого используется аппроксимация уравнений методом опорных операторов. Для решения системы алгебраических уравнений используется предобуславливапие (ILUO, ILUK, ILUTP) и различные виды итерационных методов (сопряженные градиенты и его модификации, метод минимальных невязок и др.). На втором этапе по полученным значениям смещений восстанавливается компоненты тензора напряжений и на основе критериев разрушения определяются возможные зоны разрушения. При достижении критического значения в критерии разрушения, счет прекращается. Тем самым, непосредственно процессы разрушения не моделируются.
В процессе гидроразрыва, часть жидкости проникает в разлом CD, провоцируя высокие напряжения у концов разлома. Спровоцированные напряжения в разломе могут привести к его росту. Разлом заполнен более сыпучей и проницаемой породой. Вследствие этого, а также проникновения жидкости в разлом его берега начинают движение по касательной и по нормали, что демонстрируют графики (рис.3.8), при различных углах разлома.
Как видно, разлом начинает „раздуваться", его объем при этом увеличивается и берега движутся в противоположные стороны. Возникают усилия, приложенные к берегам по нормали, которые могут привести к разры ву породы в местах локализации напряжений. Тангенциальные смепдоиия имеют различные значения на берегах, тем самым вызывая проскальзывание берегов оіішсательпо друг друге. Максимальные напряжения на на разломе достигаются вблизи его концов, что видео на рисунках (рисЗ.П).
При проведении гидравлического разрыва наииолее внжиан "ихцача, это калача определения нон разрушения. Обычно , щт моделировании процес -еов гидравлического разрыва оласш рассматриваю-! двумерные пли пеов-др-грехмерные постановки. Задача решаеіся относительно переменных: раскрытие трещины и ее длина. При етевд предполагается, что трещина развивается еамоподобоо. П такой ситуации трещина не может менять направление своего роста, что сомнительно. В более общей постановке направление возможного роста трещины может быть определено с шшошыо критериев разрушении, например критерия Кулона-Мора. Трещина может развиваться по нескольким основным сценариям, нто и демонстрируют рисунки. В зависимости от угла наклона разлома, вида приложенных нагрузок к трещине гидрорачрыва, начального распределения полей напряжении направленно возможного роете трещины будет различно. Это определяется видом возникающих напряжений и взаимодействием зон: локализации напряжений.
