Многозонный метод граничных элементов и его применение к задаче инициации трещины гидроразрыва из перфорированной обсаженной скважины Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 16, № 6, 2011

Многозонный метод граничных элементов и его применение к задаче инициации трещины гидроразрыва из перфорированной обсаженной скважины*

Д. В. Есипов, Д. С. Куранаков, В.Н. Лапин, С. Г. Чёрный Институт вычислительных технологий СО РАН, Новосибирск, Россия

e-mail: yes . еер. of f @gmail. com

Представлен многозонный метод граничных элементов решения трехмерных уравнений линейного упругого равновесия для кусочно-однородных тел (рассматривается внешняя задача упругости), заключающийся в применении к каждой из подобластей с разными упругими параметрами метода граничных элементов и в последующем объединении полученных матриц в результирующую матрицу псевдожесткости. Проведено тестирование разработанного метода на задачах, имеющих аналитическое решение. Метод применен для моделирования инициации трещины из перфорированной обсаженной и необсаженной скважин в породе. В качестве критерия инициации трещины использовано превышение максимальным возникающим растягивающим напряжением прочности породы на разрыв. В зависимости от взаимного расположения перфорации и главных напряжений в естественном залегании породы получено критическое давление жидкости, при котором инициируется трещина гидроразрыва. Установлено, что наличие обсадной колонны значительно повышает давление, необходимое для инициации трещины.

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

Введение

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

Бурение скважины производится в породе, находящейся в уже напряженном состоянии. При этом давление в скважине поддерживается равным некоторому давлению бурения Р^, необходимому для предотвращения разрушения стенок скважины. В скважину с давлением Р^ устанавливается стальная обсадная колонна. Затем производится перфорация обсадной колонны и породы, что позволяет уменьшить давление разрушения. Эффективность перфорирования зависит от распределения напряжений в естественном залегании породы, ориентации перфорации и других факторов [1-3]. Далее

*Работа выполнена при финансовой поддержке РФФИ (грант № 11-01-00475-а).

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

В работе предлагается численная модель инициации трещины гидроразрыва из перфорированной обсаженной скважины. Основной частью модели является расчет напряженно-деформированного состояния породы в окрестности обсадной колонны и перфорации, Модель также включает критерий разрушения, в качестве которого используется превышение критического напряжения на разрыв породы ас. Зоны разрушения полностью определяют начальную конфигурацию трещины гидроразрыва. Давление, при котором образуются зоны разрушения, есть давление инициации гидроразрыва.

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

Однако в случае кусочно-однородных материалов такой метод непосредственно не применим, так как при построении фундаментального решения, лежащего в основе метода, используется предположение об однородности материала в рассматриваемой области. Предложенный в настоящей работе многозонный метод граничных элементов позволяет разрешить внешнюю задачу упругости для кусочно-однородных тел, состоящих из множества однородных тел с известными параметрами упругости. Метод заключается в применении к каждой из подобластей метода граничных элементов. После аппроксимации интегральных уравнений полученные матрицы объединяются в результирующую матрицу псевдожесткости, которая связывает узловые напряжения со смещениями (аналогично матрица жесткости в методе конечных элементов связывает узловые силы со смещениями),

1. Многозонный метод граничных элементов

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

1.1. Постановка задачи упругости

Вся область задачи упругости У = У\ + V2 разбивается та две подоблаети: У1 с границей дУ\ = ¿1 и У2 с границей дУ2 = Б2, которые связаны между собой через интерфейсную границу ¿12, как показано на рис, 1, Тогда в каждой из подобластей У\ и У2 области

Рис. 1. Двухзонная задача упругости: подобласть У1 с внутренней и внешней на бесконечности границами и подобласть У2 с границ ей 52, часть кото рой Б12 является общей с VI

выполнены уравнения статического равновесия

(х)

Е

дХг

(х)=0, х Е VI ,У2, (1)

где — компоненты тензора напряжений. Производные по пространственным координатам здесь и далее обозначаются индексом, отделяемым запятой Индексы г и ] принимают значения 1, 2, 3. По совпадающим индексам производится суммирование.

Для каждой из подобластей VI и V2 выпишем граничные условия на части границы за исключением границы между областями $12:

^ (х) = Оц (x)nj (х) = Н1 (х), х Е \^12,

и(х) = а13(х)п,(х) = Н2(х), х Е ¿2\^2. (2)

Здесь Н1 и Н2 — некоторые известные функции, а п^ — компоненты единичной нормали к поверхности. Отметим, что если область У1 не является внешней, то допускаются смешанные граничные условия на 31\Б12 и Б2\Б12 (т.е. на одной части границы задаются смещения, на другой ее части — напряжения).

На интерфейсной части границы $12 выполняются условия полной сцепки, т. е. условия равновесия сил и неразрывности перемещений:

^ (хх) = иI (х2), и (х1) + и (х2) = 0, (3)

где хх Е 51, х2 Е 52 и координаты точек хх и х2 совпадают.

Так как на бесконечности напряжения отсутствуют, то для внешней подобласти У1 добавим еще условие на бесконечности [4]

щ(ю) = 0. (4)

Дополняя уравнения (1) законом Гука с параметрами упругости Еъ и1 в области У1 и Е2, щ в области У2, получим уравнения Навье для компонент смещений. Таким образом, поставлена замкнутая дифференциальная задача упругости.

1.2. Построение результирующей системы линейных уравнений

Применим для каждой из подобластей У1 и У2 обыкновенный метод граничных элементов (подробно об этом методе см., например, [5-7]), Выпишем граничные уравнения для смещений:

с.ц(х')щ(х') = / и^(х', x.)tj(x)dS — Тц(х', х)из (x)dS, х' е вк, к = 1, 2. (5)

1

1

Здесь су (х') = 2б* (х) дая внутренней 3адачи, % (х') = — 2^ (х') для внешней задачи в случае гладкой границы в точке х', Функции и^ (х', х) и Т^ (х', х) — известное решение задачи Кельвина о действии сосредоточенной силы. Выражения для и^ и Т^ можно найти, например, в [8,9], Таким образом, уравнения (5) связывают значения компонент смещений и напряжений на границах подобластей У1 и V2.

После аппроксимации граничными элементами интегральных уравнений (5) получим следующие системы линейных алгебраических уравнений:

[ и! и12 ]( = [ X! Т12 ] [ и и21 ] ( £ ) = [ X Т„ ]

и1 и12 и2 и12

(6)

Здесь в верхнем равенстве первые п1 уравнений соответствуют участку границы ¿"Дб^, следующие п12 — участку граннцы Б12. В нижнем равенстве первые п2 уравнений соответствуют участку границы 32\312, следующие п12 — Б12. В матрицах, обозначенных и1; и2, и12, сосредоточены интегралы по ядру и^, вТ^ Т2, Т12 — интегралы по ядру Т^.

Далее обычно [6] строится объединенная система линейных уравнений, которая

и1 и2 и12 Т1 Т2 Т12 форме граничных условий (3) объединенная система принимает вид

Т1 Т12 0 0 и1 и1 и12 0 0 / tl \

0 0 Т2 Т21 и12 0 0 и2 и21 ^12

0 I 0 — I и2 0 0 0 0 t2

0 0 0 0 и21 0 I 0 I V t21 /

(7)

где I — единичная матрица. Систему (7) можно переписать в более компактном виде

Т1 Т12 0 Т21

и12 0

и21 Т2

( и! \ и!2 ^12 и2

и1 0 0 и2

()

(8)

Размерность результирующей матрицы в (8) равна п1 + п2 + 2п12.

В настоящей работе предлагается метод построения результирующей системы мень-

и

(матрицы и по построению невырождены [10]), получим для каждой из подзадач матрицы псевдожесткости К

К11 К12 и1

К13 К14 _ V и12

К21 К22 " и2

К23 К24 и21

[Их И12 ] 1 [ Т1 Т12 ] [ И2 И21 Г1 [ Т2 Т21 ]

их

и12 и2 и21

^12 ^21

(9)

Далее, используя методику, аналогичную методике объединения матриц жесткости в методе конечных элементов [11,12], и учитывая граничные условия на интерфейсной границе Б12 (3), получим систему линейных алгебраических уравнений для объединенной задачи

К11 К12 0 и1

К13 К14 + К24 К23 1 и12

0 К22 К21 и2

0 ^2

(10)

Размерность результирующей матрицы в (10) равна п1 + п2 + п12, что меньше, чем в (8), Решая систему линейных алгебраических уравнений (10), найдем компоненты смещений на 81\Б12 и 5,2\5'12 объединенной задачи упругости. Затем, используя системы (6), вычислим компоненты реакций на границе Б12. Отметим, что в случае смешанных граничных условий на 81\812и 5'2\5'12 в уравнении (10) потребуется перенести неизвестные напряжения из вектора правых частей в вектор неизвестных.

Предложенный подход позволяет решать задачи со сколь угодно большим количеством подобластей и моделировать часть подобластей задачи методом конечных элементов (для этого из матрицы жесткости, используя принцип виртуальной работы, можно построить матрицу псевдожесткости и наоборот).

1.3. Разрывные граничные элементы

Вышеописанный многозонный метод граничных элементов с аппроксимацией изопара-метрическими граничными элементами (подробнее об изопараметрических элементах см, [13]) применим в случае, если границы задачи являются гладкими, т.е. во всех точках этих границ однозначно определена единичная нормаль [14], Если же граница включает в себя участки ребер или примыканий, то существуют точки, где однозначное определение единичной нормали отсутствует. Кроме того многозонный метод граничных элементов допускает только одно значение вектора напряжений в таких точках, а следовательно, напряжения аппроксимируются непрерывно вдоль границы, в том числе через ребра и примыкания. Также важно, чтобы получаемые в методе граничных элементов матрицы были квадратными, т, е, каждой точке соответствовало только одно значение напряжения.

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

к «

точек хк внутрь элемента на некоторое расстояние в локальной системе координат

элемента , £2), как это показано на рис, 2 для линейного треугольного и билинейного четырехугольного элементов. Аппроксимация геометрии элемента и функций на разрывном элементе в этом случае будет совпадать с аппроксимацией на непрерывном элементе. Тогда координаты элемента и компоненты смещений и напряжений на нем

Рис. 2. Шаблон элемента (слева) н аппроксимация функций щж ^ на нем (справа)

можно представить в следующей форме:

м

Хг (6 ) = Е Ф™ & ^ ),

т= 1

М

т (6 ,6 ) = ^ < фт &,6),

т=1 М

и (6,6 ) = Е ^ & )• (п)

т=1

Здесь М — количество опорных точек элемента, фт — интерполирующие функции.

В треугольных элементах считается, что неизвестные заданы в некоторых точках внутри элемента (х1, х2, х3) и имеют линейную аппроксимацию вдоль элемента. Геометрия элемента описывается плоским треугольником. Таким образом, геометрия, а также значения поля смещений и поля напряжений на элементе представляются посредством трех (М = 3) линейных интерполирующих функций

Ф1 (6. ,6) = (Ас — & — 6), Ф2 (6,6) = (Ас + ¿1 — 1),

Фз (6., 6) = (Ас + 6 — 1), Ас = 1(А + 1). (12)

В четырехугольных элементах считается, что неизвестные заданы в некоторых точках внутри элемента (х1, х2, х3, х4) и имеют билинейную аппроксимацию в пространстве (£1 ), т. е. линейную аппроксимацию по любому из аргументов и £2. Аналогично геометрия элемента, поля смещений и напряжений описываются четырьмя (М = 4)

билинейными интерполирующими функциями

«■ *' = 4 (! -')(! -'), *(«-е-) = -1 (х+ ')(х -') ■

*«,.6)=4 (I+>)(I+>). *«,.«2)=-4 (I -.)(I+>). (.3)

Наряду с непрерывными граничными элементами использование разрывных элементов около ребер и примыканий позволяет задавать разрывные компоненты напряжений на них в многозонном методе граничных элементов.

2. Тестовые задачи

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

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

2.1. Задача о консольной балке

Рассматривается балка квадратного сечения, нагруженная сверху давлением Р = 1 МПа. Геометрические параметры балки: длина Е = 640 мм, сторона квадрата а = 64 мм. Материал, из которого изготовлена балка, имеет следующие характеристики: модуль Юнга Е = 200 ГПа, коэффициент Пуассона V = 0.3, что соответствует стали. Левый торец балки защемлен, остальные его поверхности ссвободны от напряжений (рис. 3).

Результаты проведенных расчетов на двух последовательностях сеток представлены на рис. 4. В первой последовательности сетки состоят из 42, 168, 672, 2688 четырехугольных билинейных разрывных элементов с показателем А = 0.7 (см. верхний ряд на рис. 4; обозначим такие элементы СЗ). Во второй последовательности сетки состоят из 84, 336, 1344, 5376 треугольных линейных разрывных элементов с показателем А = 0.7 (см. нижний ряд на рис. 4; обозначим такие элементы Т).

Аналитическое решение (из теории балок, см., например, [15]) для осевой линии балки

Ра

их(г) = 777= (¿2 - 4Ьг + 6Е2)г2. (14)

Рис. 3. Задача о консольной балке

Рис. 4. Вертикальные смещения их (слева) и изменение формы балок (масштаб смещений увеличен в 10 раз; справа)

Здесь I — главный момент инерции I = —. Сравнение проводилось для точки г =

12

1 раЬА ^ „ „

в которой смещение по х равно их(ь) =--——. В табл. 1 представлены погрешно-

8 Е1

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

2.2. Задача о составной трубе

Рассмотрим две трубы круглого сечения одинаковой длины Ь = 10 м, плотно вставленных одна в другую, как это показано на рис. 5. Внутренняя труба имеет внутренний радиус а = 5 м, наружный Ь = 8 м, внешняя труба имеет внутренний радиус Ь = 8 м, наружный с = 10 м. Внутренняя труба изготовлена из материала 1 с параметрами Е\ = 25 Па, = 0.25 (соответственно коэффициенты Ламе Ах = 10 Па, = 10 Па). Внешняя труба изготовлена из материала 2 с параметрами Ех = 50 Па, = 0.25 (коэффициенты Ламе Ах = 20 Па, = 20 Па).

Таблица 1. Вертикальные смещения их (V)

Количество и тип Аналитическое Численное Относительная

граничных решение, мм решение, мм ошибка, %

элементов

42 <Э 4.8 2.45 48.95

168 4.8 3.60 25.00

672 <3 4.8 4.26 11.25

2688 4.8 4.56 5.00

84 Т 4.8 2.41 49.79

336 Т 4.8 3.55 26.04

1344 Т 4.8 4.24 11.66

5376 Т 4.8 4.54 5.41

Рис. 5. Задача о составной трубе

Рис. 6. Радиальные смещения иг в м

Внутренняя сторона внутренней трубы нагружена постоянным давлением Р = 1 Па. Внешняя сторона внешней трубы свободна от напряжений, а торцевые концы труб свободно оперты. Результаты расчетов на последовательности сеток представлены на рис. 6. Сетки состоят из 672, 2688, 10 752 разрывных билинейных элементов с показателем Л = 0.7.

Для получения точного решения воспользуемся для каждой из труб решением задачи Ламе для цилиндрической трубы [16]. Затем, приравнивая смещения наружной поверхности внутренней трубы к смещениям внутренней поверхности внешней трубы,

получим давление, с которым одна труба давит на другую:

*=р.4

где

А =

а2 Ь

В + С"

2^1 + Ах

(15)

В =

ь

2(Ь2 - а2) дх (дх + Ах)' (а2 + Ь2) + Ах а2 Ь ^ (Ь2 + с2) + А2с2

2(Ъ2 - а2) Д!(^ + Ах)

С =

2(с2 - Ь2) д2 (д2 + А2)

(16)

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

иг(с) =

Ь2 сР0 2^2 + А2 2(с2 - Ь2) д2 (д2 + А2)'

(17)

Таблица 2. Радиальные смещения иг (с)

Количество Аналитическое Численное Относительная

граничных решение, м решение, м ошибка, %

элементов

672 0.19455 0.17715 8.94

2688 0.19455 0.18994 2.36

10752 0.19455 0.19339 0.59

которая и была принята для сравнения, В табл. 2 представлены погрешности решений, полученные на рассмотренной последовательности сеток. Как и в предыдущей задаче, виден первый порядок сходимости. Получена невысокая погрешность на достаточно грубых сетках,

3. Задача инициации трещины гидроразрыва из перфорированной и обсаженной скважины

3.1. Физическая постановка задачи

Предполагается, что порода однородна и изотропна. Ее свойства характеризуются модулем Юнга Е = 20.7 ГПа и коэффициентом Пуассона V = 0.27, что соответствует осадочным породам с месторождения Вилкое (штат Техас, США), Для такой породы прочность на разрыв ас = 3.5 МПа. До бурения вертикальной скважины напряженное состояние породы определяется тензором напряжений залегания а^, главные компоненты которого

ахх = 69 МП а, ауу = 103.5 МП а, = 103.5 МП а. (18)

Данное состояние породы считается исходным и в модели ему сопоставляются нулевые смещения. Затем производится бурение скважины, после окончания которого в скважину плотно вставляется обсадная колонна. При этом стенки скважины и стенки обсадной колонны нагружены давлением бурения Р^ = 67 МПа. Далее в скважину закачивается жидкость гидроразрыва до тех пор, пока не произойдет инициация трещины. Требуется найти давление инициации и место возникновения зародышевой трещины,

3.2. Математическая постановка задачи и ее поэтапное решение

Рассматривается перфорированная обсаженная скважина, схематично изображенная на рис, 7, а. Скважина моделируется полостью в форме цилиндра диаметром И = 10 см и достаточно большой высоты Н = 10 м. Перфорация моделируется полостью в форме цилиндра длины Ь = 10 см и диаметра <1 = 2 см, затупленного на сферу. Перфорация ориентирована в горизонтальной плоскости ху под углом ¡3 к оси ж, расположена строго по центру скважины и ее ось пересекает ось скважины, В полость скважины симметрично по ее высоте плотно вставлена стальная труба высоты Н0 = 9.5 м (параметры материала: модуль Юнга Е = 200 ГПа, коэффициент Пуассона V = 0.3) с толщиной стенки И = 5 мм, имеющая в месте расположения перфорации отверстие диаметра

Для определения напряженно-деформированного состояния породы после бурения скважины необходимо решить задачу о предварительном нагружении (рис, 7, б"), В этой задаче на всей поверхности скважины и перфорации задаются следующие напряжения:

ЬА = -Рлп + а~п, (19)

где п — внешняя нормаль к поверхности. Решение (компоненты смещений) здесь рассматривается относительно исходного предсжатого состояния породы напряжениями залегания (18), Проанализируем момент после установки обсадной колонны, пока давление в скважине равно давлению бурения. Поскольку жидкость между обсадной колонной и скважиной оказывает на ту и другую одинаковое давление, то ее воздействие

а б в

Рис. 7. Математическая постановка задачи (а) и ее разбиение на две математических подзадачи (б, в)

можно заменить условием полной сцепки (3) между стенками скважины и обсадной колонны. В рассматриваемый момент они будут действовать друг на друга с одинаковым давлением Р^. Для определения изменения напряженно-деформированного состояния породы и обсадной колонны при увеличении давления жидкости необходимо решить задачу о нагружении (рис. 7, в). В этой задаче на всей поверхности стальной колонны и перфорации задаются напряжения

tA = -(Р - Ра)п, (20)

которые соответствуют изменению давления в скважине. Здесь Р — давление жидкости гидроразрыва. Все остальные поверхности задачи свободны от напряжений, так как после установки обсадной колонны напряженно-деформированное состояние не изменилось.

Затем определяется напряженно-деформированное состояние породы относительно состояния до нагружения (18). В этом случае смещения в задаче получим путем сложения смещений, полученных при решении подзадач 1, 2, а также определяемых предсжа-тием породы. Помимо найденных величин, используя фундаментальные интегральные соотношения, можно определить тензор напряжений внутри области задачи. Однако такой подход не применим для вычисления тензора напряжений на границе задачи [17]. Поэтому после определения смещений на границе, используя закон Гука в локальной системе координат, связанной с поверхностью, определяется тензор напряжений в этой системе координат. Далее находятся все его собственные значения о\ < а2 <

В качестве критерия инициации рассматривается условие превышения максимальным растягивающим напряжением а3 напряжения на разрыв породы ас. Часть поверхности задачи, где данный критерий выполняется, назовем зоной разрушения. Давление жидкости Р поднималось до тех пор, пока не был выполнен критерий инициации. Давление Pi, при котором удовлетворялся критерий инициации, назовем давлением инициации. Для расчетов была построена сетка, состоящая из 6128 граничных элементов.

Как показали результаты расчетов, в случае обсаженной скважины в отличие от необсаженной зона разрушения всегда оказывается исключительно на поверхности пер-

Р = 0о, Р{ = 70 МПа

Р = 25о, Рг = 80 МПа

Р = 50о, Рг = 102 МПа

Рис. 8. Зоны разрушения: верхний ряд — обсаженная, нижний — необсаженная скважина

Рис. 9. Зависимость давления инициации от угла Р'. 1 — для обсаженной, 2 — для необсаженной скважины

форации (рис. 8). Из рисунка следует, что для ряда конфигураций с обсадной колонной давление инициации значительно возрастает. На рис. 9 приведена зависимость давления инициации от угла направления перфорации. При малых углах как и в случае необсаженной скважины, инициация происходит на стыке скважины и перфорации. При увеличении угла зона инициации удаляется от скважины вдоль перфорации, достигая при Р ^ 50о ее конечной части.

Заключение

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

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

Список литературы

[1] Черный С.Г., Лапин В.Н., Есипов Д.В., Куранаков Д.С. Метод граничных элементов и его приложение к задаче разрушения перфорированной скважины // Тем. сборник научных статей "Краевые задачи и математическое моделирование". Новокузнецк: Новокузнецкий филиал Кемеровского гос. ун-та, 2010. Т. 1. С. 159-168.

[2] Есипов Д.В. Моделирование процесса инициации гидроразрыва пласта методом граничных элементов // Вестник КазНУ. Механика, математика, информатика. 2010. № 3(66). С. 270-277.

[31 Hossain М.М., Rahman М.К., Rahman S.S. Hydraulic fracture initiation and propagation: Roles of wellbore trajectory, perforation and stress regimes //J. Petroleum Sci. and Eng. 2000. Vol. 27, iss. 3-4. R 129-149.

[4] Купрадзе В.Д. Методы теории потенциала в теории упругости. М.: Наука, 1963. 472 с.

[5] Rizzo F.J. An integral equation approach to boundary value problems of classical elastostatics // Quarterly J. of Appl. Math. 1967. Vol. 25. P. 83-95.

[6] Бреввия К., Теллес Ж., Вроувел Л. Методы граничных элементов. М.: Мир, 1987. 524 с.

[7] Венерджи П., Баттерфилд Р. Методы граничных элементов в прикладных науках: Пер. с англ. М.: Мир, 1984. 494 с.

[8] Тимошенко С.П., Гудвер Дж. Теория упругости. М.: Наука, 1979. 560 с.

[9] Амензаде. Ю.А. Теория упругости: Учебник для университетов. М.: Высшая школа, 1971. 288 с.

[10] Becker A.A. The Boundary Element Method in Engineering. A Complete Course. N.Y.: McGraw-Hill, 1992. 348 p.

[11] Zienkiewicz O.C., taylor R.L. The Finite Element Method. 4th ed. Vol. 2. London: McGraw-Hill, 1991. 476 p.

[12] Norrie D.H., de Vries G. An Introduction to Finite Element Analysis. N.Y.: Acad. Press, 1978. 301 p.

[13] Lachat J.C., Watson J.O. Effective numerical treatment of boundary integral equations: A formulation for three-dimensional elastostatics // Intern. J. for Numerical Methods in Eng. 1976. Vol. 10. P. 991-1005.

[14] Sladek V., Sladek J. Why use double nodes in BEM? // Eng. Analysis with Boundary Elements. 1991. Vol. 8., iss. 2. P. 109-112.

[15] Ландау Л.Д., Лифшиц Е.М. Теоретическая физика. Изд. 5-е. Т. VII. Теория упругости. М.: Наука, 2007. 264 с.

[16] Седов Л.И. Механика сплошной среды. Т. 2. М.: Наука, 1970. 568 с.

[17] Sladek V., Sladek J. Improved computation of stresses using boundary element method // Appl. Math. Modelling. 1986. Vol. 10. P. 249-255.

Поступила в редакцию 21 марта 2011 г., с доработки 19 апреля 2011 г.