МЕТОД ГРАНИЧНЫХ ЭЛЕМЕНТОВ ДЛЯ ЧИСЛЕННОГО РЕШЕНИЯ ТРЕХМЕРНЫХ ЗАДАЧ МЕХАНИКИ ТРЕЩИН Текст научной статьи по специальности «Физика»

УДК 531.875.7:514.7

МЕТОД ГРАНИЧНЫХ ЭЛЕМЕНТОВ ДЛЯ ЧИСЛЕННОГО РЕШЕНИЯ ТРЕХМЕРНЫХ ЗАДАЧ МЕХАНИКИ ТРЕЩИН

А. В. Звягин 1, Н. Н. Смирнов 1 2, Д. И. Панфилов 1, А. А. Шамина 1

1 Московский государственный университет им. М. В. Ломоносова, zvsasha@rambler.ru, mech.math.msu@inbox.ru 2 Федеральный научный центр Научно-исследовательский институт системных исследований Российской академии наук, mech.math.msu@inbox.ru

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

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

BOUNDARY ELEMENT METHOD FOR NUMERICAL SIMULATIONS OF 3-D FRACTURE MECHANICS PROBLEMS

A. V. Zvyagin 1, N. N. Smirnov 1 2, D. I. Panfilov 1, A. A. Shamina 1

1 Lomonosov Moscow State University, zvsasha@rambler. ru, mech. math. msu@inbox. ru 2 System Research institute, Russian Academy of Sciences, mech.math.msu@inbox.ru

The paper suggests a modification of boundary elements numerical method, which uses displacement discontinuity method in 3D space. The advantage of this method is in essential decrease of computational elements number in 3D space because only the surface of a fracture is placed on final elements simulating the discontinuity of elastic medium. Thus the dimension of the problem is decreased. The method proved to be effective for hydraulic fractures modeling and their interaction with natural faults in host rock formations.

Keywords: 3D space, elastic medium, fracture, stress intensity factor, boundary element method, displacement discontinuity method.

Введение. Одной из актуальных задач современной механики разрушения является задача аналитических исследований концентрации напряжений в окрестности трещин в трехмерном пространстве. Увеличение проницаемости нефтенесущих пород в окрестности добывающих скважин достигается при применении технологии гидравлического разрыва. При этом жидкость гидроразрыва закачивается в скважину под большим давлением и с большим расходом, чтобы сформировать трещины [1-2]. Для повышения эффективности операции гидроразрыва необходимо проводить предсказательное моделирование формирования трещины и введения в эксплуатацию. Роль вычислительного моделирования в данной технологии трудно переоценить, поскольку трещина гидроразрыва - один из немногих технологических объектов, который не может быть подвергнут визуальному осмотру перед введением в эксплуатацию, так как находится на глубине порядка километров под землей. Существенное влияние на режим и направление распространения трещины гидроразрыва оказывает распределение напряжений в пласте, а также наличие неоднородностей материала: зон повышенной поврежденности, малых закрытых трещин, разломов [3-5]. К настоящему времени расчеты продвижения трещин гидроразрыва были основаны на моделях либо Хри-стиановича [6], либо Перкинса (PKN) [7]. Обе модели существенно полагаются на упругую

(или пороупругую) модель для описания скелета. В рамках подобных подходов получено множество аналитических автомодельных решений [8-9]. Подробнее об используемых аналитических и полуаналитических решениях изложено в соответствующих главах работы M. J. Economides and K. G. Nolte [10].

Модели продвижения трещины гидроразрыва в упругой среде в двумерной постановке с учетом возможных неоднородностей оказались весьма эффективными благодаря простоте их вычислительной реализации. Работы [3-5, 11] положили начало серии исследований, которые привели как к получению фундаментальных результатов [2], так и к созданию коммерческих симуляторов гидроразрыва [12] в средах с естественной трещиноватостью.

В настоящее время существуют хорошо развитые эффективные методы решения двумерных задач о трещинах. Одним из таких методов является метод разрывных смещений [13]. Преимуществом данного метода является возможность точного выполнения уравнений теории упругости. При этом граничные условия выполняются на дискретном множестве точек границы, которое можно сделать сколь угодно плотным. Для трехмерных задач механики твердого деформируемого тела чаще всего используются методы конечных элементов. Но их использование в механике трещин в трехмерном пространстве сталкивается с большими трудностями, поскольку построение полей напряжений и перемещений в окрестности трещин требует построения достаточно мелкой, адаптированной к геометрии трещин, сетки из конечных элементов. При наличии системы трещин сложной геометрии задача становится фактически невыполнимой.

В данной работе предлагается численный метод граничных элементов, реализующий метод разрывных смещений в трехмерном пространстве. Преимуществом данного метода является то, что на конечные элементы разбивается только поверхность трещин, моделирующая разрыв упругой среды. Это понижает размерность задачи на стадии ее решения. С точки зрения математической теории, данный подход является одной из реализаций метода разложения решения по неортогональным функциям [14-17]. После численного определения коэффициентов разложения мы имеем фактически аналитическое представление решения в виде конечного ряда внутри области. С точки зрения памяти, нам надо хранить только найденные коэффициенты разложения, позволяющие найти любые требуемые характеристики в любой точке области решения. Это существенно с точки зрения простоты практического использования полученного решения. В отличие от предыдущих результатов [15-17], при получении которых разложение осуществлялось по функциям, являющимся потенциалами простого слоя, данная работа использует новые фундаментальные решения, являющиеся потенциалами двойного слоя.

2. Построение системы функций разложения

Основой метода разложения решения по неортогональным функциям является построение системы линейно независимых решений основной системы уравнений задачи. В статической теории упругости - это уравнения равновесия.

Введем следующие обозначения: (xl, x2, x3 ) — декартовы координаты в некоторой системе координат с базисом e, e2, e3; u. (xl, x2, x3 ), (i = 1,2,3) — компоненты вектора перемещений; X,^— упругие модули Ляме; S^, (i, j = 1,2,3) — матрицы компонент тензора деформаций и тензора напряжений; для упрощения записи воспользуемся следующим обозначением частной производной df /dxk = fk ; повторяющийся индекс в любом выражении будет означать операцию свертки; V — оператор градиента функции; V2 — оператор Лапласа.

Ограничимся случаем отсутствия массовых сил. Равновесие упругой среды описывается уравнениями Ляме:

+ !) Щ л + = 0, (I = 1,2,3). (1)

Напряжения определяются обобщенным законом Гука:

с = Х8..£,, + 2ц£ , £ = -^и + и ). (2)

Л кк г* ]г 2 ],1 1,7

В результате подстановки деформаций в закон Гука (2) получим представление напряжений с помощью перемещений

с = ^Щк,к + М(и7,1 + и17) . (3)

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

функцией У2е = 0. Применим к уравнениям равновесия (1) оператор Лапласа и примем во внимание то, что дилатация является гармонической функцией. В результате приходим к известному факту, что компоненты вектора перемещений являются бигармоническими функциями:

V2 V2 щ = 0, г = 1,2,3. (4)

Воспользуемся теоремой Альманси [18] о представлении бигармонической функции в форме

V2V2u = 0, и = (+ху, V2( = 0, VУ = 0 . (5)

Рассмотрим частное решение уравнений теории упругости в форме

У У У 4 ^

щ = х — , и2 = х — , Щ = ( + X — = 0, Ау = 0. (6)

5х1 Зх2 Зх3

Подставим (6) в уравнения Ламе (1)

1 = 1 (1 + м) (хзу,п + хзу,22 + ^3,3 + (хзу,з ) з ^ + ^У2 (хзу,1) = 0; 1 = 2, (Л + !) (х3у,11 + х3у,22 + (,3 + (х3у,3 ) 3 ) 2 + (х3у,2 ) = 0

1 = 3, (Л + !) (х3уп + х3у 22 + (,3 + (х3у,3 )3) + (х3у,3) = 0

С учетом гармоничности функций ( у, получим следующие соотношения [(1 + ц)(ъъ+(Л + 3м)у] = 0, (к = 1,2,3) . Данные равенства будут выполняться тождественно, если искомые функции связаны соотношением

У = :—. (7)

А + 3м

Таким образом, уравнениям равновесия теории упругости удовлетворяет следующее поле перемещений

и = -

X + и дф

Хт

Х + 3и дх1

и2 =-

X + и дф

Хт

X + 3и 3 дх 3

, из = фз

X + и дф

Хт

X + 3и дх3

(8)

если Аф3 = 0.

Возьмем в качестве гармонической функции ф3 потенциал двойного слоя с плотностью и( % ^)

Фз(х) = Яи(£)"

д

1

, Аф = 0,

5 дп4 |Х - ^

Фз+ (% о ) - Фз- (% о ) = 4ли(% о ), если % о е 5, 5: х3 = 0, \х\< \, |х2| < й2; |х - %| ф

(9)

х-

Полученное поле перемещений

и =_А х Ф и = -А х

дх1

дФз дх0

и = Фз - А Х

дФз дх.

(10)

где А = (X + и)/(X + зи) удовлетворяет уравнениям равновесия упругой среды (1). Для поля перемещений (10) можно определить деформации:

Б11 = АХзФз,11; ^22 = А ХзФз,22; Бзз = (1 А)Фз,з АХзФз,зз;

1,1з'

Б = - А х3ф312 ; 2б13 = (1 - А)ф - 2А х3ф 1

2Б2з = (1 - А)Фз,2 - 2А ХзФз,2з; е = % = (1 - А)фз,з-Подстановка (11) в закон Гука (2) определяет напряжения:

(11)

^ = А1 (1 - А)фз,з - АХзФз,11; = А1 (1 - А)фз,з - А ХзФз,22;

^зз _ (л , А Л л. ™ . а12

2 и

(1+ А1 )(1-А)фз,з -Ахзфз

зз'

2 и

-А х3 ф

12

(12)

сг1з (1 - А) Ст2з (1 - А)

— _ -Фз,1 - А ХзФз,1з ; ТГ1 = —Фз,2 - А ХзФз,2з •

2и 2 ^ 3,3>13 2и 2 В формулах (12) использовано обозначение А = 2и) .

Отметим, что аналогичным образом можно построить еще два частных решения уравнений упругости с полями перемещений:

и1 =ф

X + и дф

Хт

X + зи 3 дх^ 2

и = -

X + и дф X + и дф

_1_ V _'-1 11 —__1_ V _•

■ х

и

X + и дф

Хл

1 X + зи 3 дХ1

и2 =Ф2

X + зи дХ2 X + и дф

, из =-

Хт

Хл

X + зи дХ:

из

X + зи дХ3 X + и дф

Х

X + зи дХ3

1

Функции ф (Х, Х, Х3 ), (к = 1,2), как в рассмотренном случае, являются потенциалами двойного слоя (9). Для них можно определить соответствующие поля деформаций и напряжений. В результате для площадки 5 с нормалью, направленной вдоль оси Х , мы имеем три линейно независимых решения уравнений упругости. Отметим, что в случае выбора плотности потенциала двойного слоя и(ГГ) в виде многочлена, функции (9) вычисляются аналитически, поскольку сводятся к вычислению интегралов вида

гт е

ЯгГ ГГ .

У х- %

(13)

5о |х-%|з

Например, в случае постоянной плотности и прямоугольной площадки 5 : Г < \, |Г < К2, аналитическое выражение интеграла (13) будет иметь следующий вид

Х,

2 х,

:[ (-Х + (Ф, -Ф 2 -Ф5 +Ф б) + Sign (Х + К) (Фз -Ф 4 -Ф7 +ф8)]^

где

Ф =

агб

Ф =

Ф =

Ф =

агБ

агБ

агБ

sign(-х + к)(х - к х )2 + (х -к )2 + х32 - хз - ( х - к )2 + * (к2 - х2 ) | х3 |

^^ х2 * | х^ |

-sign(-х + К)(Х - К )>/( х -К )2 + (х2 -К2 )2 + хз + х32 + (х1 - К )2 + ¿'(К2 - х2) |х3

-К2 + х2 + * х3

у 2 2 , з, у

^ -sign(х + К)(Х + К)>/(х + К)2 + (х -К)2 + х32 -х2 -(х + К)2 + КК -х2)|х31^

^^ Х2 * | х^ |

¿7£п(х1 + К)(Х + К )\/(Х + К)2 + (х -К2)2 + Х3 + Х2 + (х1 + К)2 + - Х) |хз

Ф 5 = агБ

Ф =

Ф =

Ф =

агБ

агБ

агБ

-К + Х + * Х

sign(-х + К)(-х + К )>/(х -К )2 + (х + К2 )2 + х2 + х2 + (х - К )2 + КК + х) |х3

I / ^Хо

у 2 2 I Я у

^ -sign(-х + К)(-х + К )>/(х - К )2 + (х + К )2 + х2 - х2 - (х - К )2 + + Х) Х1 ^

у 2 2 , з, у

( , 1 \ /- »4 /Г . . ^ ""I \ 2 . 2 . 2 , 1 \ 2. ■, 1 \ I Л

К + х + г'|х3|

sign( х + К)(Х + К )^/(Х; + К )2 + (х + К2 )2 + х2 + х2 + (х + К )2 + г'(К2 + Х) | х3

-V

-х2 +1 хз

у

-sign(х + К)(х + К )>/(х + К )2+(х + К )2 + х3 - х2 - (х + К )2 + г'(К2+Х) |х3

К+х + * х3

Для построения численной схемы обозначим поля перемещений и напряжений, соответствующие каждому из построенных решений, следующим образом:

Ц?к\о*к\ (и = 1,2,3), (т = 1,2,- • (14)

В выражениях (14) верхний индекс к = 1,2,з означает номер решения. Например, номер к = з соответствует решению (8), (12). Индекс т - это номер граничного элемента.

Численный метод решения краевой задачи. Рассмотрим типичную задачу механики трещин. В пространстве глобальной системы координат ( X, У, Z ) с базисом ( зр э2, э3)

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

(Хт, У,, )— координаты центра граничного элемента с номером т , (е^, е%, ет ) — локальный базис данного элемента. Введем в рассмотрение ортогональную матрицу ат, поз-

17

воляющую выразить локальные векторы в глобальном базисе

(15)

а7 э7 .

Для каждого граничного элемента с номером п в его локальном базисе еП можно считать

~ п

заданными три компоненты вектора напряжений на площадке с нормалью е :

С31 = Ь1П ,С32 = Ь2 ,С33 = Ь3 . (16)

Рис. 1. Система трещин в глобальной системе координат На поверхности одной из трещин выделен граничный элемент со своей локальной

системой координат ( х1, х2, х3 ) в базисе е, е2, е3

Рассмотрим следующее поле перемещений и напряжений

ит = Дтит(1) + о2тит (2) + Вътит(Ъ)\

ут _ ГА т т (1) 7-у т т (2) , Р) т _т (3)

Лу = и1 + и2 + и3 с1] ■

(17)

Эти поля являются решениями уравнений теории упругости в локальном базисе е™,е™,еът .

Найдем координаты центра элемента с номером п в данном базисе с использованием матрицы перехода(15)

хт" =( Хп — Хт) ат +(Уп — Ут) а12 +(2 — ¿т) а^,

х'т" =( Хп — Хт ) а™1 +(Уп — Ут ) ат22 +(¿п — ) а£, (18)

хтп =(Х — Х )ат +(У — У )ат +(2 — 2 )ат■

3 \ п т; 31 \ п т / 32 \ п т / 33

Подстановка координат (18) в формулы (17) позволяет найти компоненты тензора напряжений от поля (17) в точке, соответствующей центру элемента с номером п

^т ( тп тп тп\ _ тлт т(1)/ тп тп л-тп\

(, л ^ , х3 ) — С (, х ^ , х3 ) +

7 7 (19)

гл т т(2)/ тп тп л-тп\ т~у т т(3)/ тп тп л„тп\

+ ^2 (х1 ,х2 ,х3 ) + ^3 (х1 ,х2 ,х3 ).

Полученные компоненты тензора (19) в базисе в™,в™,в™ позволяют найти их значения в

базисе в",в2 ,в

где

<У,

тп _ т-ч т тп(1)

У

Пт

7

т тп(2)

У

т

+А (

т тп(3)

т

+ Б3 Оу

_тп(1) _т (1)

7 = 7

У РЧ

_тп(2) _т(2)

7 . = 7

У РЧ

7У

тп(3) _ т(3)

сг

РЧ

(атап)(атап),

\ Р8 иД дг ]Г )•>

(атап V атап),

\ Р8 вД дг ]Г }•> ()().

\ Р^ иД дг уг)

(20)

Выражения (20) принято называть коэффициентами влияния напряжений, которые вызывают соответствующее единичное решение для элемента с номером т в центре элемента с номером п . Причем напряжения (20) соответствуют базису в ", в2 , в3". Это позволяет найти вклад в граничные условия в напряжениях в форме девяти матриц размерности (Ж X N), где N — общее количество граничных элементов

()(<х), ()(атк),

()(«таг),

Л(к) тп тп(1) = 73к = 7т(1) РЧ

В(к) тп тп ( 2 ) = 73к = 7 т(2) РЧ

С(к) тп тп (3) = 7зк = 7 т(3) РЧ

к = 1,2,3.

(21)

Коэффициенты О™, Б™, Б™ (их количество 3М) определяются из граничных условий. Суммируя с неопределенными коэффициентами вклады (21), которые вносит каждый элемент в центр элемента с номером п , и используя граничные условия, получим

У Л(1)БТ + У В(1)БТ + У С(1)= Ь",

/ 1 тп 1 / * тп 2 / * тп 3 1 '

т=1 т=1 т=1

МММ

У Л(2)+ У В(2)+ У С(2)= Ьп,

тп 1 тп 2 тп 3 2

т=1 т=1 т=1

N N N

У Л(3)БГ + У В(3)+ У С(3)= Ьп,

тп 1 тп 2 тп 3 3

т=1

т=1

т=1

(22)

>

Введем общий вектор из коэффициентов и вектор правой части (22) в следующей форме Если сформировать глобальную матрицу коэффициентов влияния в виде

~А(1) тп В(1) тп С (!) тп

м = А2) тп В(2) тп С (2) тп

Л(3) тп В(3) тп С (3) тп

то система уравнений (22) перепишется в векторной форме

(23)

МБ' =Ь'

(25)

где ^,Ъ' — транспонированные векторы (23).

Таким образом, задача сводится к решению линейной системы уравнений (25). Если эта задача решена, то определение перемещений и напряжений в любой точке глобальной системы координат сводится к следующей последовательности действий:

В цикле по переменной т = 1...^ определяем перемещения щ(Х,У,2) и напряжения с (х,у,г).

1. Находим координаты данной точки в локальной системе координат с базисом

„т Л т Л т /г-ч

^ , е2 , е3 с использованием матрицы перехода (5)

т аи а 21 а31

т Х2 = а12 а22 а32

т V Хз У _ а13 а2з а33

'X — ХтЛ

т

У - У

т

7 — 7

v 7 7 т у

(26)

2. Вычислим величину перемещений и напряжений

ит, Ът ,(1 = 1,2,3; , = 1,2,3)

с помощью формул (17), (26).

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

т Т тт т

и, = и а

т т т т

Ъ..а а. = с .

V Щ ,Р РЧ

4. Суммируем в цикле полученные вклады

N

N

иг(Х,у,г) = Xит, с,(Х,у,г) = ХЪ,. (27)

т=1 т=1

Таким образом, формулы (27) завершают определение перемещений и напряжений в произвольной точке области решения.

3. Результаты расчетов и тестирования

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

Моделирование дискообразной трещины. Рассмотрим осесимметричную трещину в виде диска [19-21], которая находится под действием внутреннего давления.

В цилиндрической системе координат г ,ф, г (трещине соответствует диск г = 0,0 < г < Я ) данной задаче соответствуют граничные условия: 2 = 0, 0 < г < Я, с„=— р, с = 0.

? ? 22 г ? г2

Основной характеристикой линейной механики разрушения является коэффициент интенсивности напряжений К = с^2я(г — Я) . В аналитическом решении для

круглой трещины его величина равна

К =■

2 ЯГС22 (г^ г

{ V я 2 — 1

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

Таблица 1

Радиус трещины Я (м) 20 25

Численное значение К1 № ( м12) 0,0063 0,0074

Теоретическое значение К /№ ( м12) 0,0051 0,0056

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

(а) (б)

Рис. 2. Форма цилиндрической трещины:

а) соответствует раскрытию трещины радиуса Я = 20 с прямоугольными граничными элементами

К х к, ^ = 1К =1,

б) два сечения Х = 0 и у = 0 верхней поверхности трещины радиуса Я = 5 , = 0,35, Л = 0,35

Плоская трещина под действием неравномерного внутреннего нагружения. В качестве еще одного сравнения была выбрана плоская трещина, которая находится под действием линейно возрастающего внутреннего давления Р = Р( х) (рис. 3).

Рис. 3. Плоская трещина длины 21 под действием линейно распределенного давления

Р(х), Р( А) = 0, Р(В) = р

Теоретические значения коэффициентов интенсивности напряжений в данном случае соответственно равны: в точке В К = 3ру[л1 /4 ; в точке А К = р4л1 /4 (рис. 3). Для сравнения были проведены расчеты трехмерной трещины в виде плоского прямоугольника размерами Ь х Ь = 8,4 х 2,4 . Граничные элементы определялись параметрами

X у

Нх = 0,1, Ну = 0,1. Коэффициенты интенсивности определялись для сечения у = 0, соответствующего середине трещины. Результаты сравнения полученного численного решения с точным аналитическим решением данной задачи приведены в табл. 2.

Таблица 2

K / Л ( м12 ) в точке A K /( (м|2 ) в точке B

Расчет 0,00053 0,00149

Теория 0,00048 0,00139

Сечения верхнего берега (z = 0+ ) трещины приведены на рис. 4.

Берег трещины в сечении у= ■ 0.1 (середина Берег трещины в сечении х= - 0.1 (середина

1 -0.? 0 0.? 1 -4 -2 О 2 4

координата X координата Y

(а) (б)

Рис. 4. Сечения верхнего берега (Z = 0+ ) плоской трещины, находящейся под действием нарастающего по оси х внутреннего давления:

а) сечение y = 0; б) сечение X = 0

Представление трехмерной поверхности трещины кривыми ее пересечения с плоскостями y = const приведено на рис. 5.

Рис.4б показывает выполнение симметрии задачи в сечении x = const. Линейно распределенное давление внутри трещины приводит к асимметричному раскрытию трещины (рис. 4а, рис. 5). Расчетные значения коэффициента интенсивности (табл. 2) достаточно хорошо соответствуют их теоретическим значениям.

Рис. 5. Поверхность верхнего берега трещины представлена кривыми пересечения верхнего ( 2 = 0+ ^

берега трещины с плоскостями У = СОПБ1

4. Решение трехмерных задач взаимодействия трещин

Для демонстрации пространственных возможностей программы приведены иллюстрации расчетов двух задач (рис. 6, 7).

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

Раскрытие пассивной трещины

Рис. 6. Две круговых плоских трещины Я = 0,1, Я = 0,05 расположены под углом 30°

на расстоянии 0,01:

(а) относительное расположение; (б) раскрытие пассивной трещины

В первой задаче активная трещина радиуса Я = 0,1 находится под действием внутреннего

давления р/№ = 0.005 (рис. 6). Под действием возникающих напряжений пассивная (разрез)

трещина радиуса Я = 0,051 раскрывается. На рис. 6б показан график разности нормальных

перемещений верхнего и нижнего берега вспомогательной трещины в сечении Х = 0.

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

Рис. 7. Круглая в плане трещина радиуса Я = 0.05 изогнута в виде цилиндрической поверхности радиуса Я = 0.025:

а) общий вид;

Ь) перемещение берегов трещины по нормали к серединной поверхности в сечении х = 0

Внутри трещины действует давление Я/2№ = 0,001. Рис. 7Ь - график зависимости нормальных перемещений нижнего и верхнего берега трещины от дуговой координаты в сечении Х = 0. Следует отметить наличие морщин на внешней поверхности трещины.

Рис. 8. (а) плоскость сечения трещины; (Ь) график изменения второго инварианта девиатора тензора напряжений в данном сечении

Для иллюстрации распределения напряжений в окрестности трещины было определено значение второго инварианта девиатора тензора напряжений (рис. 8Ь) в зависимости от координат (X, г) в сечении у = 0 (рис. 8а) и аналогичный график от координат (у, г) (рис.

9Ь) в сечении X = 0 (рис. 9а).

Сечен1ккр\-спсп1111П11ндр11ческсп1треи1инып-юскоапъюх = 0 Поле значений второго инварианта девиатора напряжений

Рис. 9. (а) плоскость сечения трещины;

(Ь) график зависимости второго инварианта девиатора тензора напряжений в данном сечении

Приведенные расчеты фиксируют, что образование морщин приводит к существенному снижению уровня напряжений (рис. 8Ь), поскольку в данном сечении у = 0 трещина фактически закрыта (это сечение соответствует значению дуги 5 = 0 (рис. 7Ь). В сечении X = 0 (рис. 9а) края трещины раскрыты, что приводит к бесконечным напряжениям (рис. 9Ь). Наличие морщин отчетливо проявляется в немонотонном поведении напряжений на круглом контуре трещины в данном сечении (рис. 9Ь).

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

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

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

Работа выполнена при поддержке гранта РФФИ 16-29-15076.

Литература

1. Киселев А. Б., Захаров П. П. К численному моделированию динамики необратимого деформирования и разрушения нефтеносного пласта // Матем. моделирование. 2013. Т. 25. № 3. С. 62-74.

2. Smirnov N. N., Kiselev A. B., Nikitin V. F., Smirnova M. N., Tyurenkova V. V. Underground Hydraulic Fracturing Technology Computer Simulations // Proc. The IACGE International Symposium on Geotechnical and Earthquake Engineering (IACGE-2016). Beijing, China. October 11-13, 2016. Beijing. 2016. P. 194-202.

3. Акулич А. В., Звягин А. В. Взаимодействие трещины гидроразрыва с естественной трещиной // Изв. Рос. Академии наук. Механика жидкости и газа. 2008. № 3. С. 104-112.

4. Акулич А. В., Звягин А. В. Численное моделирование распространения трещины гидроразрыва // Вестн. Моск. ун -та. Сер. 1. Математика. Механика. 2008. № 1. С. 43-49.

5. Богданов А. И., Звягин А. В., Тьерсилен М. Взаимное влияние системы трещин на коэффициент интенсивности напряжений // Вестн. Моск. ун-та. Сер. 1. Математика. Механика. 2004. № 6. С. 44-49.

6. Khristianovich S. A., Zheltov Y. P. Formation of vertial fractures by means of highly viscous liquid // Proc. 4th World Petrol. Congr. Rome, 1955. V. 2. P. 579-586.

7. Perkins T. K., Kern L. R. Widths of hydraulic fractures // J Petrol Technol. 1961. V. 13. № 9. P. 937-949.

8. Smirnov N. N., Tagirova V. P. Self-Similar Solutions of the Problem of Formation of a Hydraulic Fracture // Porous Medium. Fluid Dynamics. 2007. Vol. 42. № 1. P. 60-70.

9. Smirnov N. N., Tagirova V. P. Problem of Propagation of a Gas Fracture // Porous Medium. Fluid Dynamics. 2008. Vol. 43. № 3. P. 402-417.

10. Economides M. J., Nolte K. G. Reservoir Stimulation. Third Edition. Chichester : Win-heim : NY : Brisbane : Singapore :Toronto : John Wiley & Sons Ltd, 2000. 750 p.

11. Chuprakov D. S., Akulich A. V., Siebrits E., Thiercelin M. Hydraulic Fracture Propagation in a Naturally Fractured Reservoir // SPE128715 Presented at the SPE Oil and Gas India Conference and Exhibition held in Mumbai, India. January 20-22, 2010. Mumbai, 2010.

12. Weng X., Kresse O., Chuprakov D., Cohen C.-E., Prioul R., Ganguly U. Applying complex fracture model and integrated work flow in unconventional reservoirs // Journal of Petroleum Science and Engineering. 2014. 124. Р. 468-483.

13. Крауч С., Старфилд А. Методы граничных элементов в механике твердого тела. М. : Мир, 1987. 328 с.

14. Алексидзе М. А. Решение граничных задач методом разложения по неортогональным функциям. М. : Наука, 1978. 352 c.

15. Shou K. J. A high order three-dimensional displacement discontinuity method with application to bonded half-space problems : dis. ... Ph.D. Minneapolis : University of Minnesota, 1993.

16. Shou K. J., Siebrits E., Crouch S. L. A high order displacement discontinuity method for three-dimensional elastostatic problems // Int J Rock Mech Min Sci. 1997. № 34 (2). Р. 317322.

17. Wu K. Numerical Modeling of Complex Hydraulic Fracture Development in Unconventional Reservoirs : dis. ... Ph.D. Austin : University of Texas at Austin. 2014.

18. Almansi E. Sull'integrazione dell'equazione differenzigle v u ~ 0 // Ann Mat. 1899. V. III. № 2.

19. Уфлянд Я. С. Интегральные преобразования в задачах теории упругости. Л. : Наука, 1967. 402 с.

20. Гольдштейн Р. В. Плоская трещина произвольного разрыва в упругой среде // Изв. АН СССР. Механика твердого тела. 1979. № 3. С. 111-126.

21. Kassir M. K., Sih G. G. External elliptical crack in elastic solid // The international Journal of Fracture Mechanics : Wolters-Noordhoff Publishing Groningen. 1968. Vol. 4. № 4. P. 347-356.