CC BY
31
УДК 519.64
Д. С. Куранаков, Д. В. Есипов, В. Н. Лапин, С. Г. Черный
Институт вычислительных технологий СО РАН пр. Акад. Лаврентьева, 6, Новосибирск, 630090, Россия
kuranakov@ict.sbras.ru
ТРЕХМЕРНЫЙ ДУАЛЬНЫЙ МЕТОД ГРАНИЧНЫХ ЭЛЕМЕНТОВ РЕШЕНИЯ ЗАДАЧ УПРУГОСТИ С ТРЕЩИНАМИ *
Разработан трехмерный дуальный метод граничных элементов решения задач линейной упругости для тел с регулярной границей и с трещинами. Классический метод граничных элементов заключается в численном решении граничного интегрального уравнения смещений на регулярной границе упругого тела. В случае тела с трещинами уравнение смещений вырождается на трещиноватой границе. В дуальном методе граничных элементов предлагается на трещиноватой части границы решать граничное интегральное уравнение напряжений. Интегралы в этом уравнении имеют высокий порядок сингулярности и понимаются в смысле главного значения Коши и Адамара. Для их вычисления используется метод выделения сингулярности (разложение в ряд подынтегрального выражения и вычисление специальным способом главных значений интегралов от сингулярных членов разложения). Дуальный метод граничных элементов верифицирован на задаче о внешнем нагружении бесконечного материала с плоской круглой трещиной. Метод показывает высокую точность в вычислении полей смещений и напряжений в материале, а также коэффициентов интенсивности напряжений на фронте трещины.
Ключевые слова: линейная упругость, метод граничных элементов, трехмерное моделирование, трещина, коэффициенты интенсивности напряжений.
Введение
Одним из методов численного решения задач упругости является метод граничных элементов (МГЭ). Его главным отличием от метода конечных элементов и метода конечных разностей является отсутствие необходимости аппроксимации всей области упругого тела. МГЭ требует аппроксимации только границы области.
Суть классического МГЭ заключается в переходе от дифференциальных уравнений упругости к граничному интегральному уравнению смещений. При решении этих уравнений отыскиваются все неизвестные на границе области. Для нахождения неизвестных внутри области применяется фундаментальное интегральное уравнение для внутренней точки. Однако классический МГЭ не применим для задач с трещинами, потому что граничное интегральное уравнение смещений вырождается на берегах трещины.
Для преодоления указанной трудности в настоящей работе разработана модификация классического МГЭ - дуальный метод граничных элементов (ДМГЭ). Суть ДМГЭ состоит в замене граничного интегрального уравнения смещений на граничное интегральное уравнение напряжений на одном из берегов трещины. Это уравнение не включает в себя смещения
* Исследование выполнено при финансовой поддержке Новосибирского технологического центра компании «БсЫитЪе^ег».
Куранаков Д. С., Есипов Д. В., Лапин В. Н., Черный С. Г. Трехмерный дуальный метод граничных элементов решения задач упругости с трещинами // Вестн. Новосиб. гос. ун-та. Серия: Информационные технологии. 2015. Т. 13, вып. 1. С. 74-90.
ISSN 1818-7900. Вестник НГУ. Серия: Информационные технологии. 2015. Том 13, выпуск 1 © Д. С. Куранаков, Д. В. Есипов, В. Н. Лапин, С. Г. Черный, 2015
на трещиноватой границе, но позволяет найти разрывы смещений на трещине. Интегральное уравнение напряжений содержит особенности более высокого порядка, чем интегральное уравнение смещений, а интегралы, входящие в него, понимаются в смысле главного значения Адамара. Ввиду этого требуется создание численного метода для их эффективного вычисления.
Авторами данной работы классический МГЭ применялся для решения задачи зарождения трещин гидроразрыва в окрестности перфорированной скважины [1; 2]. Двумерный ДМГЭ применялся в задаче распространения трещин гидроразрыва от круглой полости [3]. Целью данной работы является разработка трехмерного ДМГЭ, который в дальнейшем может быть применен для решения трехмерных задач распространения трещин гидроразрыва.
Постановка задачи упругости
для бесконечного тела с полостью и трещиной
Задача упругости решается в бесконечной области V , которая имеет границу
дУ = Б = Б' + Б" + 5 ++ 5-,
как изображено на рис. 1.
Б"
Рис. 1. Внешняя задача упругости для тела с полостью и трещиной
В каждой точке х е V выполняются уравнения упругого равновесия
да, (х)
у = дх} у
= ау,,(х) = (1)
Здесь а, - компоненты тензора напряжений. Индексы 7, у принимают значения 1, 2, 3.
Индекс, отделённый запятой, здесь и далее обозначает дифференцирование по соответствующей координате. По повторяющемуся индексу производится суммирование. К уравнениям (1) добавим закон Гука для изотропного однородного материала
а у = Л5у£кк + 2 № у . (2)
Здесь £у = 0.5("7у + ",7) - компоненты тензора деформаций, и7 - компоненты вектора смещений, 5, - символ Кронекера, Л и и - коэффициенты Ламе, выражаемые через модуль
Юнга Е и коэффициент Пуассона у как
Еу Е
Л =-, и =-. (3)
(1 + у)(1 - 2у) 2(1 + у)
Добавляя к дифференциальным уравнениям (1) краевые условия на полости
(,(х )| = '*(х),
I5 * (4)
"7(Х )|Б" = "7 (Х К
на берегах трещины
а также условие на бесконечности
(х)Г + =
Г * (5)
*г (х)|Г - = АД х),
и, (х)|Г „ = 0, (6)
получаем замкнутую дифференциальную задачу. Здесь = ОуПу - компоненты вектора напряжений, пу - компоненты вектора нормали к границе области. Функции, отмеченные * -
это известные распределения смещений и напряжений.
Отметим, что уравнения упругого равновесия (1) с учётом закона Гука (2) и связи тензора деформаций со смещениями являются уравнениями эллиптического типа для компонент смещений.
Классический Метод Граничных Элементов
Одним из методов решения дифференциальной задачи (1)-(6) является метод граничных элементов (МГЭ) [4]. Его основная идея заключается в переходе от дифференциальной постановки (1)-(6) к граничному интегральному уравнению относительно неизвестных функций на границе. В классическом МГЭ рассматривается тело без трещин (Г+ = Г- = 0). Тогда граничное интегральное уравнение смещений (ГИУС) для точки коллокации у на регулярной границе Г * = Г* + Ги запишется в виде
Су (У )иу (у) = | и у (у, х)*; (х)<яТ (х) - у.р. | Ту (у, Х)и у (х)йГ (х) . (7)
Г* Г*
Здесь у - точка коллокации, х - точка интегрирования. Функции и у (у, х) и Ту (у, х) -
фундаментальное решение Кельвина [5] задачи о действии сосредоточенной в точке у силы
/*(х) = 5(х - у)е, . Коэффициент Су (у) в общем случае зависит от геометрии поверхности в
точке у . В случае внешней задачи и гладкой поверхности он равен -5у / 2 . Заметим, что для
внешней задачи интегральное уравнение (7) автоматически удовлетворяет граничному условию (6) на бесконечности [6]. Здесь и далее сингулярные интегралы у.р. и у.р.Ь рассматриваются в смысле главного значения Коши и Адамара соответственно.
Неприменимость классического МГЭ для тел с трещинами
Рассмотрим тело, имеющее регулярную границу Г* и трещину Г + + Г- . Запишем ГИУС для точки коллокации у + е Г + на одном из берегов трещины [7]
Су (у +) • и (у +) + Су (у -) • и (у -) =
= { иу (у+,х)*, (хМГ(х) - { Ту (у+,х)и, (хМГ(х) +
Г* Г*
+ | иу (у+,х)*, (хМГ(х) - у.р. | Ту (у+,х)и, (х)^Г(х) + (8)
Г + Г+
+1 и у (у -, х)*, (х ^Г (х) - у.р. | Ту (у -, х )и, (х ^Г (х).
Г- Г-
Здесь и далее индексами + и - обозначены точки верхнего Г + и нижнего Г- берегов трещины.
На трещиноватой границе Г1 введем разрыв смещений
Ли, (х1) = и, (х+) - и, (х-). (9)
Предположим, что напряжения ', на берегах трещины Б + и Б равны по модулю и противоположны по знаку
(х +) + (х -) = 0. (1°)
В задачах нагружения трещины давлением жидкости выражение (1°) не является существенными ограничением, так как давление на противоположные берега трещины одинаково, а сила трения жидкости мала по сравнению с давлением. Заметим, что функции ядер на верхнем и нижнем берегах трещины обладают следующими свойствами
и 1} (у, х +) = и 1} (у, х -),
- (11) Ту(y,х+) = -Ту(y,х ).
Тогда с учётом (9)-(11) уравнение смещений (8) запишется в виде
Су (у +) • и (у +) + Су (у -) • и (у -) = { иу (у +, х )', (х № (х) -
С + г + (12)
-{ Ту(у +,х)и,(х(х) -у.р. \Ту(у +,х)Аи,(х^Б(х).
Б* Б+
В точках у- на противоположном берегу трещины Б- ГИУС совпадает с ГИУС (12) на
границе Б +. Таким образом, ГИУС вырождается на трещиноватой границе и классический МГЭ не применим для решения задач с трещинами. Для преодоления этой трудности в [8] предложена модификация классического МГЭ - дуальный МГЭ (ДМГЭ).
Дуальный метод граничных элементов
В ДМГЭ на регулярной границе решается граничное интегральное уравнение смещений, а на трещиноватой границе - граничное интегральное уравнение напряжений.
Граничное интегральное уравнение смещений
Для точек у на регулярной границе Б * запишем ГИУС, аналогичное ГИУС (7) в классическом МГЭ
Су (у)"у (у) = | и у (у,х)'у (х^Б(х) - У.р. | Ту (у,х)"у (х^Б(х) - | Ту (у,х)Аи, (х^Б(х). (13)
Б Б Б+
Граничное интегральное уравнение напряжений
Для решения проблемы с вырождающимся ГИУС (12) на трещиноватой границе в ДМГЭ предлагается на одном из берегов трещины вместо ГИУС (12) решать граничное интегральное уравнение напряжений (ГИУН)
(у+) = [ Ц(у +,х)'у (х^Б(х) - М у (у +,х)"у (х^Б(х) - Ь.у.р. { М у (у +,х)Аи,(х^Б(х). (14)
Б* Б* Б+
Здесь Ьу(у +, х) = Вку(у +, х)пк(у +) и М у(у+, х) = Бку(у +, х)пк(у+), а функции Вкг] и Бк1] получаются из ядер и у и Ту путём дифференцирования по соответствующей координате и применения закона Гука (2) (полные выражения для функций и у, Ту, Вку , Бку можно найти, например, в [7]). Это уравнение не содержит компонент смещений и1 на трещиноватой границе, но позволяет найти неизвестные компоненты разрывов смещений Аи, на границе. Главной сложностью ДМГЭ является создание эффективной и точной процедуры вычисления главного значения Адамара интеграла в уравнении (14).
Дискретизация границы
и получение системы линейных алгебраических уравнений (СЛАУ)
Продемонстрируем численный метод решения ГИУН (14) на примере пространства с трещинами (Я = S+ + S—). Вся граница области S аппроксимируется граничными элементами Яе, как показано на рис. 2 слева,
Я - I^. (15)
е=1
Рис. 2. Разбиение границы Я на разрывные квадратичные граничные элементы Яе (слева) и параметризация (^1,^2) элемента с Ма = 9 (справа)
Каждый граничный элемент Яе параметризуется локальными координатами (^1,^2), как показано на рис. 2 справа. Параметризация проводится таким образом, что любой криволинейный четырехугольный элемент трансформируется в стандартный квадрат с вершинами (1;1), (1;—1),(—1;—1),(—1;1), а криволинейный треугольный элемент - в стандартный треугольник с вершинами (0,0) ,(0,1),(1,0). Компоненты радиус-векторов, разрывов смещений, напряжений в произвольной точке элемента представляются в виде
= 1 /А^Ш^Ъ). (16)
а=1
Здесь ха - узловые точки элемента, фа(4\,42) - интерполирующие функции элемента, Ыа -количество узлов и интерполирующих функций в элементе. Уравнение (14) с учётом (15) и (16) запишется в виде
NeNa( e ) ( ^
(17)
tj(y") = Z Z ~AueaiiMff(y"¿¿Ш^Й)JШгЖМг e=1 a=1 ^ £ &
Здесь J(¿,¿2) - якобиан перехода в локальную систему координат элемента. Разрывы смещений Auf1 в a -й узловой точке e -го элемента выносятся за знак интегрирования, так как не зависят от переменных интегрирования и ¿2. Заметим, что интегралы в формуле (17) зависят только от геометрии границы области и не зависят от граничных условий.
Записав уравнение (17) в узлах yea, получим СЛАУ относительно неизвестных функций Au i:
MAu = -t. (18)
Здесь Au и t - векторы, составленные из компонент разрывов смещений и напряжений во всех узловых точках; M - матрица, составленная из значений интегралов в уравнении (17).
В случае тела с регулярной границей и трещинами (Б = Б шется в виде
Тп - 0.51
42
м
21
м
22
Г и > Аи
и
11
21
0 ' -I
Г1 ^
V1 у
Б + + Б ) система (18) запи-
(19)
где и и Т - подматрицы, составленные из значений интегралов по элементам в ГИУС (13), Ь и М - подматрицы из ГИУН (14).
Граничные элементы и аппроксимирующие функции
Так как ГИУН (14) требует гладкости поверхности в точках коллокации у на трещине
Б1, в ДМГЭ необходимо использовать разрывные элементы, все узлы которых находятся внутри элемента, как показано на рис. 2. В настоящей работе использовались разрывные линейные и квадратичные элементы [7], а также специальные элементы для фронта трещины [9], определяемые набором интерполирующих функций
-?)(. -1)(. - р)
Фу(ЛхЛг) =
=
2Л2(д -1)(д - р) -(# -Щх +?)(. -1)(. - р) т]2(д -1)( д - р)
Ф(Р Р ) = £(£ +Л)(. -1)(. - Р) 2л (д -1)(д - р)
Ф^ъЬ) = Ф(6,&) = Ф6(Ь,4г) =
^7(^1,^2) =
Ф8(#1,#2) =
Ф9(#1,#2) =
£(£ +у)(. - д)(. - р) 2^2( д -1)( р -1) ' -£(£ +?)(. - д)(. -1) 2Л2(д - р)(р -1) ' -?)(£ + ?)(. - д)(. -1) ^2( д - р)( р -1) -£(£ -?)(. - д)(. -1) 2^2( д - р)( р -1) '
-^)(. - д)(. - р) 2^2( д -1)( р -1) ' -(#1 - + у)(. - д)(. - р)
^2( д -1)( р -1) :
(20)
где . = 1 + , р = л/ 1 + ^ , д = ^11 , ^ - параметр, регулирующий положение узлов в элементе. В интерполирующие функции (20) входит выражение . = ^ 1 + , которое хорошо
аппроксимирует асимптотику разрывов смещений Аи на фронте трещины, что положительно влияет на точность вычисления коэффициентов интенсивности напряжений.
Вычисление главного значения Адамара сингулярного интеграла
Главная сложность ДМГЭ заключается в построении алгоритма вычисления главного значения Адамара сингулярного интеграла по граничному элементу Бе, содержащему точку коллокации у
I у (у) = Ь.У.р. | Ку (у, х)йБ(х), (21)
где подынтегральное выражение выглядит следующим образом
Б
Ку (у, х) = Му (у, X)фа(х)7(х) = Як1] (у, х)пк (у)фа (х)7(х). (22)
Для вычисления интеграла I у (21) используем метод выделения сингулярности [10]. Этот
метод заключается в разложении в ряд подынтегрального выражения (22) в окрестности особой точки у и в вычислении специальным способом интегралов от первых членов этого разложения. Таким образом, интеграл I у (21) расписывается в виде суммы трёх интегралов
— / 3 (23)
где
I у (у) = I—Чу) +1—2( у) +1—3( у)
Г
I—— Чу) = | Ку (у, х) —
Я
I-—2( у) = У.р. |
ку; 2(в)
к—2(0) Кг—3(0)^
р
р
& (х),
р
2*Мв) ку2{в) dS(х) = Нш I I ———рйрйв,
0 ре(0)
р
ЦЧу) = Ь.р.у. I dS(х) =1хш
р
£^0
г 2„мв) к^{в) н у) ^
| | -Jр-Lрdрdв +
0 ре(в) р 8
(24)
(25)
(26)
Здесь р и в - локальные полярные координаты элемента с центром в точке коллокации у , изображенные на рис. 3. Функции К—2(в) и К—3(в) - первые члены разложения в ряд по переменной р подынтегрального выражения К у (22), рЕ(в) - радиус £ -окрестности в локальной полярной системе координат элемента, р*(в) - расстояние до границы элемента. Заметим, что £ -окрестность в локальной системе координат не является окружностью и рЕ(в) зависит от полярного угла в . Коэффициент Н(у) в (26) характеризует асимптотику
в точке у бесконечной части интеграла I—Чу) [12].
Я
Я
и в локальной системе координат (Ъ,Ъ ) (справа)
Вычисление слабо сингулярной части интеграла
Для вычисления слабо сингулярного интеграла Цу1 = О (1 / р) (24) в настоящей работе был
применен метод трансформации области интегрирования [11]. В этом методе элемент разбивается на треугольники с вершинами в точке сингулярности О с локальными координатами (^1,^2), как показано на рис. 4, а. Заменой переменных (Ъ,Ъ2) ^ (£1,£2) треугольник ОАВ трансформируется в стандартный треугольник (рис. 4, б).
ъ -1 ъ-1 Ъ + 1 ъ -1
(£ - 1) - (#2 - 1) (£ - 1) + (#2 + 1)
(27)
Стандартный треугольник заменой переменных (£1,£2) ^ (ух, у2) трансформируется в квадрат [-1;1] х [-1;1] таким образом, что одно из его ребер ( у2 = -1, у1 е [-1;1]) соответствует особой точке О (рис. 4, с). Преобразование квадрата в треугольник (у1,у2) ^ (£1,£2) выглядит следующим образом:
£(У1, У2) = 4(1 + УХ)(1 + У2), Í2(v1,У2) = 1(1 - У1)(1 + У2).
(28)
Тогда обратное преобразование из треугольника в квадрат (£1,£2) ^ (у1, у2) запишется в виде
^^ & + &2 * 0 Ь1 + &2
[[-1;1] & +&2 = 0 (29)
[2(& +&2) -1 +&2 * 0 [-1 +&2 = 0
Якобиан трансформации треугольника в квадрат = (1 + у2)/8 = О(р) аннулирует сингулярность 1/р и позволяет вычислять интегралы от слабо сингулярных функций (К/р)
К 1 1
[ = [ [ —
оАв Р -1 -1 р
(30)
Рис. 4. Разбиение элемента на треугольники (а) трансформация треугольника ОАВ в стандартный треугольник (Ь) и в квадрат с ребром О (с)
На рис. 5 изображен характерный вид слабо сингулярной функции, возникающей в ДМГЭ, и квадратурная сетка для ее численного интегрирования.
Рис. 5. Распределение подынтегрального выражения К33 - К332 / р2 - К333 / р3 по поверхности граничного элемента и квадратурная сетка его для численного интегрирования
Вычисление главных значений сильно сингулярных и гиперсингулярных интегралов
Для вычисления оставшихся двух сильно сингулярного I-2 = О (1/ р2) (25) и гиперсингулярного I- = О(1/ р3) (26) интегралов необходимо вычислить члены К~2 и К- разложения по р подынтегрального выражения в (21), а также разложить в ряд по £ радиус рЕ(в).
Для этого запишем первые члены разложения радиус-вектора г = х - у по переменной р
г(в) = а(в)р + Ь(в)р2 + О(р3), (31)
где коэффициенты а(в) и Ь(в) выражаются через угол в и производные радиус-вектора х по локальным координатам £ и £2 в точке коллокации у
а(в) = х,£ ео8(в) + х,£ 8ш(в), (32)
Ь(в) = 2(х,££ 0О82(в) + 2 х,£1£281п(в)с08(в) + х£2£^1П2(в)). (33)
Тогда длина вектора г =| г | запишется в виде
г(в) = ра(в) (1 + рВ(в) + О(р2)), (34)
где
а(в) =|а(в)|, (35)
в(в) = <аввр. (36)
Запишем уравнение (34) для радиус-вектора £ -окрестности и выразим из него радиус р£(в)
р (в) = £-^--£2Вв! + О(£3). (37)
£ а(в) а 2(в) (37)
Проинтегрируем аналитически сингулярные интегралы по переменной р
2л
у = Нш I К——2(в)(1п(р*(в)) — 1п(р£(в)))в :
2л
2л
2л
(38)
| Кг;—2(в)1п(а(в)р*(в))dв + liш01n£ | К~2(в^в = | К—2(в)1п(а(в)р*(в))в,
IГ3 = Нш
7 £—>0
Г 2л
(I
V 0
I К— 3(в)
1
1
рн(в) р*(в)
К у)
2л
= Нш
£— 0
2л
IК 3(в
V 0
ав + вв"
£ р*(в)
К у)
(39)
I К— 3(в) Vв(в)—'
0
Заметим, что коэффициенты разложения К~2(в) и К~3(в) зависят только от точки кол-локации у и переменной в и не зависят от переменной р, что позволяет провести интегрирование по р аналитически. Значения контурных интегралов (38) и (39) подсчитываются численным интегрированием.
Разложение в ряд подынтегрального выражения
Ядро Яку в (22) имеет вид [7]
где
Яку (у, х) = лЬ) ((у, х) + у, х) + Я§?( у, х) + Я(4>( у, х))-)—,
Яй)(y, х) = 3гтПт ((1 — 2^уГк + у(3гкГ,у + §}кг,,) — 5ГгуГк ), Я8)( y, х) = 3у(п ,Г угк + пуГ г,к X Я§)( y, х) =(1 — 2у)(3пкК Г у + пуАк + п ёук X Я§)( у, х) = —(1 — 4у)пк8у,
г = Г
,1 г
(40)
(41)
(42)
(43)
(44)
(45)
Для разложения в ряд по р ядра Яку (40) необходимо найти члены разложения от его составных частей
1
(3рв(в) + О (р-1))
где
г 3( у, х) р3^ (в)
г/у,х) = А (в) + рdR1 (в) + О(р2),
а, (в)
А (в) =
а(в)
dRI (в) = — А (в) В(в). а(в)
Вектор нормали в точке интегрирования х представляется в виде
где
п, (х) = п1 (у) + рdNI(в) + О (р2), dNi(в) = пъ ео8(в) + пц 8т(в).
(46)
(47)
(48)
(49)
(50)
где
м_^ (в)
Б-уу (в) = ГТТ-^Г (67)
4^(1- у) а3(в)
2 м Б\у (в)-3В(в) 51, (в)
Б-2 (в) =-£--ку ' У ' ку '. (68)
к,уУ ' 4^(1- у) а (в)
Здесь
Б°к„ (в) = Бку1)(в) + Бку (в) + Б^в) + Бк,у4)(в), (69)
(56)
Разложим в ряд функции 5$ (41), Бк2) (42), Бк33) (43), 5$ (44). Первая функция 5^ (41) запишется как
Бк;/)( у, х) = Бк0у1}(в) + рБЦ71)(в) + О(р2), (52)
Бй1)(в) = (аи (в)пи (у ))г1у (в), (53)
Б^в) = 2]у (в) [аи (в)пи (у)] + 2% (в) (в)п„ (у) + ат {в)<Шт (в)]. (54)
Здесь
2% (в) = 3[(1-2у)Ак (в)ду +у(Ау (в)8л + А, (в)у )-5Ак (в)Л, (в)Ау (в)], (55)
21у. (в) = 3[(1 - 2у)ёЯк (в)5у + уЩ вк + Ж, (в)5}к) --5^ (в) Л (в) Ау (в) + Ак (в)Щ (в) Ау (в) + Ак (в) Л (вЩ (в))]. Вторая функция Б^ (42) запишется как
Б$( у, х) = Б°ку\в) + рБ1Р(в) + О (р2), (57)
где
Б°ку\в) = 3у[п, (у) Ау (в) Ак (в) + пу (у)А (в) Ак (в)], (58)
БЦ72)(в) = 3у[п, (у)(Ау (в)ёЯк (в) + Ак (вЩ (в)) + пу (у)(А (вЖ (в) + Ак (в)Ж, (в)) ■
+ёНг (в) Ау (в) Ак (в) + ёНу (в) Аг. (в) Ак (в)]. Третья функция 5^3 (43) равна
Б§)( у, х) = Бкуз)(в) + рБ1(у3)(в) + О (р2), (60)
где
Бкуз)(в) = (1-2у)[3пк (у) А (в) Ау (в) + пу (у )дк + п, (у )5]к ], (61)
Б^в) = (1 - 2у)[3пк (у)( А (в^Яу (в) + Ау (в)ёЯг (в)) + +3ёНк (в) А (в) Ау (в) + (в)5,к + ёнг(в)5]к ]. Четвертая функция Б^у (44) запишется как
Б$( у, х) = Б0') +рБу + О (р2), (63)
где
БкЦ4)(в) = -(1-4у)пкду, (64)
Б1у4)(в) = -(1-4у)ёНкду. (65)
Тогда с учетом (46)-(65) разложение ядра Бку (40) выглядит как
Бку (у, х) = Б- (в)р-3 + 5-2 (в)р-2 + О (р-!), (66)
где
(59)
(62)
Slv (в) = S^V)+Sf(e)+Sf(e)+Sf(e). (70)
Для разложения в ряд подынтегрального выражения (22) необходимо вычислить члены разложения функций фа(x) и J(x).
Интерполирующие функции фа в точке интегрирования x запишутся в виде
фа(х) = фа (у) +^Ф(в) + O(p\ (71)
где
d Ф(в) = фа
Л\ cos(e) + фаЛ2 sin(e). (72)
Якобиан J в точке интегрирования х равен
J (х) = J (у) + pdJ (в) + O (p2), (73)
где
dJ(в) = J,£ cos(e) + J^ sln(e). (74)
Тогда разложение в ряд подынтегрального выражения Ку (22) в запишется в виде
К (у, х) = К-3(в)р"3 + К- 2(в)р~2 + Op1), (75)
где
ку 3(в) = nk (у) S-3 (в)фа( у) J (у). (76)
К- 2 (в) = nk (у )[Ski2 (в)фа(уУ(у) + S-3 (в)(фа(у )dJ(e) + dФ(в)J(y))]. (77)
Итак, мы получили два первых члена К- и К- разложения подынтегрального выражения K j (22) в ряд по локальной координате p .
Метод экстраполяции смещений
для вычисления коэффициентов интенсивности напряжений
Для вычисления коэффициентов интенсивности напряжений (КИНов) в точке O фронта трещины, как показано на рис. 6, использовались интерполяционные формулы [7]
К? = —V < , (78)
7 4(1 - у )\ 2l b
-? = E
4(1 -y2)\2l
Kn = Ы ~A<, (79)
K°in —Auf, (80)
111 4(1 + k)V2Z " V 7
где AuP - разрыв смещений в точке P трещины на расстоянии l от точки фронта ° , векторы b, n и t - базис локальной системы координат на фронте трещины. Уравнения (78), (79), (80) известны как одноточечные интерполяционные формулы. Можно заметить, что их точность напрямую зависит от точности вычисления полей смещений в окрестности фронта трещины. Эти формулы применимы, если расстояние l достаточно мало по сравнению с характерным размером трещины. Если же l велико, значения КИНов получаются заниженными. В этом случае необходимо использовать экстраполяцию значений КИНов KP1 and KP2, посчитанных по двум точкам P1 и P2 на расстояниях l1 и l2 соответственно, на точку ° фронта трещины (рис. 6).
K° = KP2 +12( ^-KPi), (81)
l2 l1
Уравнение (81) известно как двухточечная интерполяционная формула.
Рис. 6. Вычисление КИНов методом экстраполяции смещений в окрестности фронта трещины
Верификация ДМГЭ
Для верификации ДМГЭ рассмотрена следующая задача. В бесконечном материале, растягиваемом в направлении у напряжениями а, имеется плоская круглая трещина радиуса Я, повернутая вокруг оси О2 на угол а , как показано на рис. 7.
О
а
а >
/ ^^ч
/ \
X
\Я
Рис. 7. Задача о круглой трещине, наклоненной на угол а относительно растягивающих напряжений а
Ширина трещины Ж в этой задаче описывается соотношением [13]
Ж(г) = асо82а8(1 у ^Я2-г2, (82)
лЕ
а КИНы на фронте трещины равны ([14; 15])
2Я
К1 = 2асо8 а—, (83)
4 „ 1Я
Кп =-азтасозасозв.—, (84)
2- у Мл
— 4(1- у) . Я
Кш =-а81пасо8а81пв./—, (85)
2-у V л
где в - угловая координата, характеризующая положение точки на фронте трещины.
Для верификации метода построена серия сеток, характеризующихся разбиением в окружном направлении на Ыв элементов и в радиальном - на Ыг элементов. Брались значения физических параметров Я = 1т, р = 1МРа, Е = 2°GPa , у = 0.2 , а = 0 .
На рис. 8 показана сеточная сходимость численного решения задачи к точному.
Рис. 8. Распределение раскрытия трещины Ж по всему радиусу г (слева) и в окрестности фронта трещины (справа): 1 - Ыг = 8 (пунктир); 2 - 16 (пунктир-точка); 3 - 32 (точка),
сплошная - точное решение. Щ = 32
Видно, что ДМГЭ вычисляет раскрытие трещины Ж с достаточно высокой точностью. Максимум погрешности достигается на фронте трещины. Порядок аппроксимации составляет ~1/2.
Посмотрим, как погрешность вычисления раскрытия трещины Ж влияет на погрешность в вычислении КИНов.
На рис. 9, 10 приведены значения К1 во всех узлах сетки вдоль одного радиуса, посчитанные по формуле (78) для К1 (рис. 9 - линейные элементы, рис. 10 - квадратичные).
Рис. 9. КИН К1 : точное решение (83) (пунктир); формула (78) для точного решения смещений (82) (сплошная линия) и для численных решений ДМГЭ с линейными элементами: Ыф =32; О - Ыг =8,
□ - 16 , д - 32
1.4
1.3
1.2
1.1
0.9
0.8
0 7 С_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I
' 0 0.2 0.4 0.6 0.8 1
г,ш
Рис. 10. КИН К1 : точное решение (83) (пунктир); формула (78) для точного решения смещений (82) (сплошная линия) и для численных решений ДМГЭ с квадратичными элементами: Nф = 16; О - Nг = 4, □ - 8, д -
16
1.4
1.3
1.2
1.1
0.9
0.8
0.7
0.2
0.4 0.6
г,ш
0.8
0
Рис. 11. КИН К1 : точное решение (83) (пунктир); формула (78) для точного решения смещений (82) (сплошная линия) и для численных решений ДМГЭ с квадратичными элементами и специальными элементами (20) на фронте трещины): Nф = 16 ; О - Nг = 4 , □ - 8, д - 16
Видно, что на расстоянии нескольких элементов от фронта К1 вычисляется достаточно точно, а на фронте трещины погрешности велики. Замена линейных элементов на квадратичные не дает никакого улучшения точности (см. рис. 10).
Рис. 10. Зависимость КИНов от положения точки фронта плоской круглой трещины: точное решение (сплошная); К[ (О ), Кп (□), Кш (д)
Рассмотрим сетки, аналогичные сеткам на рис. 10, но с прифронтовыми элементами специального типа (20). Результаты расчетов K¡ приведены на рис. 11.
С использованием специальных элементов КИН K¡ вычислен с погрешностью не более 2 % на самой грубой сетке (Nr = 4 ).
о
Рассмотрим задачу о круглой трещине, наклоненной на угол а = 45 относительно растягивающих напряжений. На рис. 12 показаны распределения всех трех КИНов по фронту трещины.
КИНы вычислены по формуле (81). Погрешность вычислений составляет не более 2 %.
Заключение
Разработан трёхмерный дуальный метод граничных элементов для решения внешних и врутренних задач упругости для тел с трещинами. Для вычисления гиперсингулярных интегралов применен метод выделения особенности. Для повышения точности вычисления КИНов на фронте трещины применены элементы специального типа, учитывающие особенность поведения решения на фронте трещины. ДМГЭ верифицирован на задаче о нагруже-нии плоской круглой трещины и показывает высокую точность в вычислении полей смещений и напряжений в материале, а также коэффициентов интенсивности напряжений на фронте трещины. Разработанный метод может быть применен в задачах распространения трещин гидроразрыва.
Список литературы
1. Есипов Д. В., Куранаков Д. С., Лапин В. Н., Черный С. Г. Многозонный метод граничных элементов и его применение к задаче инициации трещины гидроразрыва из перфорированной обсаженной скважины // Вычислительные технологии. 2011. Т. 16. № 6. С. 13-26.
2. Alekseenko O. P., Potapenko D. I., Cherny S. G., Esipov D. V., Kuranakov D. S., Lapin V. N. 3D Modeling of fracture initiation from perforated non-cemented wellbore // SPE J. 2013. Vol. 18, No. 3. P. 589-600.
3. Алексеенко О. П., Есипов Д. В., Куранаков Д. С., Лапин В. Н., Черный С. Г. Двумерная пошаговая модель распространения трещины гидроразрыва // Вестн. Новосиб. гос. ун-та. Серия: Математика, механика, информатика. 2011. Т. 11, № 3. С. 36-59.
4. Rizzo F. J. An Integral Equation Approach to Boundary Value Problems of Classical Elastostatics // Quart. J. of Applied Mathematics. 1967. Vol. 25. P. 83-95.
5. W. Thomson (Lord Kelvin), On the equations of equilibrium of an elastic solid // Cambr. Dubl. Math. J. 1848. Vol. 3. P. 87-89.
6. Купрадзе В. Д. Методы теории потенциала в теории упругости. 1963. М.: Наука. 472 с.
7. Aliabadi M. H. The Boundary Element Method: Vol. 2. Applications in Solids and Structures. 2002. John Wiley and Sons Ltd., 598 p.
8. Mi Y., Aliabadi M. H. Dual boundary element method for three-dimensional fracture mechanics analysis. Engineering. Analysis. 1992. Vol. 10. No. 2. p. 161-171.
9. Cisilino A. P., Aliabadi M. H. Three-dimensional BEM analysis for fatigue crack growth in welded components // Int. J. for Pressure Vessel and Piping. 1997. Vol. 70. P. 135-144.
10. GuiggianiM., Krishnasamy G., Rudolphi T. J., Rizzo F. J. A general algorithm for numerical solution of hypersingular equations // J. Appl. Mech. 1990. Vol. 57. P. 906-915.
11. Lachat J. C., Watson J. O. Effective numerical treatment of boundary integral equations: A formulation for three-dimensional elastostatics // Int. J. for Numerical Methods in Engineering. 1976. Vol. 10. P. 991-1005.
12. Hadamard J. Lectures on Cauchy's problem in linear partial differential equations. New Haven: Yale University Press, 1923. 316 p.
13. Sneddon, I. N., Elliott, H. A. The Opening of a Griffith Crack Under Internal Pressure // Quart. Appl. Math. 1946. Vol. 4, P. 262.
14. Murakami Y.(Editor-in-Chief) Stress Intensity Factors Handbook. Pergamon Press, 1987.
15. Tada H., Paris P., Irwin G. The Stress Analysis of Cracks Handbook. ASME Press, NY 2000. 3rd ed.
Материал поступил в редколлегию 15.12.2014
D. S. Kuranakov, D. V. Esipov, V. N. Lapin, S. G. Cherny
Institute of Computational Technologies SB RAS 6 Lavrentiev Ave., Novosibirsk, 630090, Russian Federation
kuranakov@ict.sbras.ru
THREE-DIMENSIONAL DUAL BOUNDARY ELEMENT METHOD FOR ELASTICITY PROBLEMS WITH CRACKS
Three-dimensional dual boundary element method for elasticity problems with regular boundary and with cracks is developed. Conventional boundary element method consists in numerical solution of displacement boundary equation on regular boundary of the elastic body. In case of the body with cracks displacement boundary equation degenerates on the crack surface. Therefore in dual boundary element method it is suggested to solve traction boundary integral equation on the crack surface. Integrals in this equation have the singularity of high order and are considered as Cauchy and Hadamard principal values. Substraction of singularity is used for their evaluation. This method consists in derivation of Lourent series of the functions under integral sign and special method of integration of the first members of the series. Dual boundary element method is verified by the problem of penny-shaped fracture in unbounded material. The method shows high accuracy in calculation of stress and displacement fields, as well as stress intensity factors at the front of the fracture.
Keywords: linear elasticity, boundary element method, three-dimensional modeling, crack, stress intensity factors.
