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

УДК 532.591, 531.5.031

ВЗАИМНОЕ ВЛИЯНИЕ ДИСКООБРАЗНЫХ ТРЕЩИН В ТРЕХМЕРНОМ УПРУГОМ ПРОСТРАНСТВЕ

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

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

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

The mutual influence of three-dimensional disc-shaped cracks located in parallel planes of an elastic medium is studied. The medium is under tensile stress directed perpendicular to the crack planes. The cracks are modeled by mathematical cuts of the continuous medium with the possibility of a strong rupture of the displacement field on the banks of the cut. The discontinuous displacement method is used to solve the problem numerically.

Key words: three-dimensional space, elastic medium, crack, stress intensity factor, boundary-element method, discontinuous displacement method.

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

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

1 Звягин Александр Васильевич — доктор физ.-мат. наук, проф. каф. газовой и волновой динамики мех.-мат. ф-та МГУ; гл. науч. сотр. ИМАШ им. A.A. Влагонравова РАН, e-mail: zvsasha®rambler.ru.

2 Панфилов Дмитрий Игоревич — канд. физ.-мат. наук, инж. каф. газовой и волновой динамики мех.-мат. ф-та МГУ, e-mail: panfilovdmQmail.ru.

1 Шамина Анастасия Александровна — асп. каф. газовой и волновой динамики мех.-мат. ф-та МГУ, e-mail: anashamina90® mail .ru.

поля перемещений и напряжений, возникающие в упругой среде. Важную роль в механике трещин играет величина коэффициентов интенсивности напряжений (КИН) [3] на краю трещины. Часто в качестве критерия начала роста трещины в определенном направлении принимается достижение коэффициентом интенсивности (или комбинацией коэффициентов интенсивности напряжений) некоторого критического значения, характерного для данного конкретного упругого материала. Поэтому одной из составляющих задачи является проверка метода на возможность определения КИН и сравнение его численных значений с имеющимися аналитическими результатами.

"°'5 0 "0,5 -1 -0,5 оТГГТГо "°'5

Рис. 1. Исследуемые типы расположения двух трещин в пространстве: а центры основной и дополнительной (малой) трещин лежат на одной оси: б малая трещина расположена над основной трещиной:

в малая трещина не затеняет основную трещину

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

3. Построение системы функций разложения. Основой метода разложения решения по неортогональным функциям является построение системы линейно независимых решений системы уравнений задачи [4|. В статической теории упругости это уравнения равновесия. Введем следующие обозначения: (х1,х2,хз) — декартовы координаты в некоторой системе координат с базисом в1, е2, е3; и^(х1, х2, х3), г = 1, 2, 3, — компоненты вектора перемещений; А, ц — упругие модули Ламе; ai,j, г,] = 1, 2, 3, — матрицы компонент тензора деформаций и тензора напряжений; для упрощения записи воспользуемся следующим обозначением частной производной: ^^ = /д.; повторяющийся индекс в любом выражении будет означать операцию свертки; V — оператор градиента функции; V2 — оператор Лапласа.

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

(А + ц) пк,ы + цuijj = 0, г = 1,2,3. (1)

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

а3г = <Ц,7 7.7.- + 2/Х^г, £ = \ + Щ^)

(2)

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

aji = А5jiUk,k + Ц (uj,i + и^-).

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

У2УЧ

0,

г = 1, 2, 3.

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

V2V2ui = 0, и = + х3^, VV = 0, V2ф = 0.

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

дф дф дф ^ ,

, и2 = хзт—, и3 = <рз+жзт—, V Уз = 0, У^ф = 0. (3)

0Ж1 дж2 джз

Подставив (3) в уравнения Ламе (1), будем иметь

г = 1, (А + (жзф,и + жзф,22 + Уз,з + (жзф,з),з)д + ^2(жзфд) = 0; г = 2, (А + (жзф,и + жзф,22 + Уз,з + (жзф,з),з);2 + ^2(жзф,2) = 0; г = 3, (А + (жз^,ц + жз^,22 + Уз,з + (жзф,з),з);з + ^2(жзф,з) = 0.

С учетом гармоничности функций Уз, ф получим следующие соотношения:

[(А + ^)уз,з + (А + 3^)ф,зЦ = 0, к = 1,2,3. Они будут выполняться тождественно, если искомые функции связаны соотношением

А + ^

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

А + ^ <9уз А + ^ <9уз А + ^ <9уз

Щ = --——— ж3 , и2 = - ж3 , мз = Уз ~ х , » х3 , (4)

А + 3^ дж1 А + 3^ дж2 А + 3^ джз

если V2у3 = 0.

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

¥»(*> ^У^ = 0, (5)

в любой точке поверхности 5, где 5 : жз =0, |ж1| ^ |ж2| ^ Л,2, обладающего свойством у+ ({о) — ^з (Со) = 4тф({0), если {о € 5, причем уз ~ на бесконечности. Полученное поле перемещений

щ = -Лж3^, и2 = -Лжз^-, = Уз - Лж3^, (6)

дж1 дж2 джз

где Л = удовлетворяет уравнениям равновесия упругой среды (1). Для поля перемещений (6)

можно определить деформации

£11 = —Лжзуз,11, ^22 = —ЛжзУз,22, £зз = (1 — Л) уз,з — ЛжзУз,зз,

£12 = —ЛжзУз,12, 2£1з = (1 — Л) Уз,1 — 2ЛжзУз,1з, (7)

2£2з = (1 — Л) Уз,2 — 2ЛжзУз,2з, е = = (1 — Л) Уз,з.

Подстановка (7) в закон Гука (2) позволяет определить напряжения 711 722

— = Л1 (1 - Л) у3;3 - Лжзуздь = Л1 (1 - Л) у3)3 - Лж3уз,22,

азз (712

= (1 + Л1) (1 -Л)у3;3 -ЛжзУз.зз, = -Лж3уз,12, (8)

(713 (1 - Л) Д 723 (1 - Л) Д

7Г~ = —«—¥>3,1 - ЛжзУздз, — =---Уз,2 - Лж3уз,23,

2^ 2 2^ 2

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

«1 =

Л + ^

-ж3-

«1

Л + 3^ дж1 Л + ^

-хз-

Л + 3^ дж1

Л + ^

и-2 = ~ х , о , «з

Л + 3^ дж2

Л + ^

г/,2 = <Р2- \ , о , из

Л + 3^ дж2

Л + /л д<р\ _ Л + 3^ дж3' Л + ^

-хз-

Л + 3^ дж3

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

^3(ж) =

йо

ст си 41 ; 42

|ж-£|;

^£1 ^£2-

(9)

Например, в случае постоянной плотности, при т = п = 0, и прямоугольной площадки £о |£1| < ^1, |£21 < Л-2 аналитическое выражение интеграла (9) будет иметь следующий вид:

-/ е Жз(Ж1 - £1) Жз(Ж2 - £2) , Дж1 - £1)(Ж2 - £2)

/(Ж1,Ж2,Ж3,£1,£2) = —;—:-—--;—;-— + аг^.ап 1

Г(Г + Ж2 - £2) Г(Г + Ж1 - £1)

а?з(а?1 ~£1)(ж2 ~Ы(х'1 + г2) г(г2ж| + (ж! - £1)2(ж2 - £2)2)'

ж3г

где г(ж1,ж2,ж3,£ь£2) = л/{х\ - б)2 + (ж2 - £2)2 +

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

и

т(к)

а,

т(к) Ч '

г,; = 1,2,3, т = 1,2,Ж.

Здесь верхний индекс к = 1, 2, 3 означает номер решения. Например, к = 3 соответствует решет

4. Численный метод решения краевой задачи. Рассмотрим типичную задачу механики трещин. В пространстве глобальной системы координат (X, У, Я) с базисом (пъП2, П3) расположено несколько трещин, которые моделируются поверхностями разрыва перемещений (рис. 2). Будем для определенности ставить задачу в напряжениях. Это означает, что на берегу трещины задан вектор напряжений как функция точек поверхности. Пусть (Хт, Ут, Ят) — координаты центра граничного элемента с номером т; (ет, ет, ет) — локальный базис данного элемента. Введем в рассмотрение ортогональную матрицу а™, позволяющую выразить локальные векторы в глобальном базисе:

-0,5

ег •

(10)

п

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

аз1 = Ьи, аз2 = ^ азз = Щ •

(11)

Рис. 2. Система трещин в глобальной системе координат (X, У, Я). На поверхности одной из трещин выделен граничный элемент со своей локальной системой координат (ж1, Ж2, жз) в базисе ех, е2, е3

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

ит_ит(1) + + ^™ит(з)

гт _ г-)т7т(1) _1_ Г>т,-_т(2) I т~\т—т(з) ( )

= 7у' + ^2 + ^з 7у ,

где — неопределенные коэффициенты.

Эти поля являются решениями уравнений теории упругости в локальном базисе (еЗ^еЦ, еЦ). Найдем координаты центра элемента с номером п в данном базисе с использованием матрицы перехода (10):

ж]?т = (Хп — Хт) ап + — а12 + — "оз

ж™га = (Хп — Хт) а21 + — а22 + — ^т) а2з' (13)

= (Хп — Хт) а31 + — а32 + — азз-

Подстановка координат (12) в формулы (11) позволяет найти компоненты тензора напряжений

п

Ет (™тга ™тга ™тга) _ пт7т(1) (жтп ™тга ™тга) + пт7т(2) (жтп ™тга ™тга) +

¿7 (ж1 ,ж2 ,жз )= 7 ¿7 (ж1 ,ж2 ,жз ) + 7 ¿7 (ж1 , ж2 , жз ) +

(з) (жГ, ж11",Ж1") . (14)

С помощью полученных компонент тензора (13) в базисе (еЗ^еЦ, е—) можно отыскать их значения в базисе (е", е", е"):

где

тга _ г>т,-г™га(1) I г>т7™га(2) | г>т7™га(з)

7у = 7у + + ^з ,

т"(1) _ ™(1) т„га ) ^тп(2) _ т(2)

,т „га

7у = 7Р<? "г ^ , 7у = 7Р<? ^"дг"г

7тга(з) = 7т(з) а„ ) / т а„

7= 7Р9 1ар«аг^ I адг"г

(15)

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

■ А (к) _ тга(1) _ т(1) / тга ) ( т п )

Атга = 7зк = 7рд арз"з^ I адг%г/ ,

в!" = 7зтга(2) = 7ртд(2) (а^) ( а1 а"г) , (16)

с(к) = 7т"(з) = 7т(з) ( ата " ) ( ата " ) стга = 7зк = 7рд ^ ар5аз^ ^адг акг/ '

Коэффициенты (их количество 3Ж) определяются из граничных условий. Сум-

мируя с неопределенными коэффициентами вклады (16), которые вносит каждый элемент в центр п

( N N N

Е ^т + е ^т + ^ с« ^т = ^,

т= 1 т= 1 т= 1

NN N

Е

т= 1 т= 1 т= 1

+ е ^т + ^ с(2" ^т = (17)

=1 =1

N N

А&я т + Е ^т ^ с(з"^т = ь?, п = 1,2, N.

NN N

Е

ч т= 1 т= 1 т= 1

Введем общий вектор из коэффициентов:

Б= DN, ^ 1, Ь, ^ N) (18)

и вектор правой части (17):

ь = ь, ^, , Ь2, ь, ^, , Ь2, ь, ^).

Если сформировать глобальную матрицу коэффициентов влияния в виде

М

Гд1 B1 C1

-^mn ^ mn mn a2 B2 C2

mn mn mn

A3 в3 С3

mn mn mn

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

MD* = b

(19)

(20)

(21)

где О*, Ь* — транспонированные векторы (18), (19). Таким образом, задача сводится к решению линейной системы уравнений (21). Если эта задача решена, то определение перемещений и напряжений в любой точке глобальной системы координат (X, У, 2) сводится к такой последовательности действий.

1. В цикле по переменной т = 1,..., N определяем перемещения щ (X, У, 2) и напряжения ац (Х,У,2).

2. Находим координаты данной точки в локальной системе координат с базисом б™, б™, б™, воспользовавшись матрицей перехода (16):

aii &21 031 ai2 &22 &32 й!3 a23 ^33

х — Xm

Y Ym Z — Zm

3. Вычислим величину перемещений и напряжений Ц™, Х22, г = 1, 2, 3 ] = 1, 2, 3, с помощью формул (7), (11), (20).

т

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

m ij

rm.m ) i aij ,

- a - a ' -'iju'iqu'jp

a,

pq-

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

Ui (X, Y, Z) =

N

E

m=1

U m ^ t 1

aij (X,Y,Z) =

N

E£m-

m=1

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

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

1. Осесимметричная трещина в виде диска [4-8], которая находится под действием внутреннего давления. В цилиндрической системе координат r, z (трещине соответствует диск z = 0, 0 ^ r ^ R) в данной задаче имеем граничные условия: z = 0, 0 ^ r ^ R, azz = —p arz = 0. Основной характеристикой линейной механики разрушения является КИН К \ = lim 1тгг^/2п(г — R).

r^R

В аналитическом решении для круглой трещины его величина равна

R

2

Ki =

razz(r) dr

VTTRJ Л/ R2 — г2 ' 0

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

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

2. Плоская трещина, которая находится под действием линейно возрастающих) внутренних) давления Р(ж) = 0.5Р(1 + ж/1) —1 ^ ж ^ I. Теоретические значения КИН в данном случае соответственно равны К \ = в точке х = I и К\ = 'рл/тй/А в точке х = —I. Были проведены расчеты трехмерной трещины в виде плоского прямоугольника размера х = 8.4 х 2.4. Граничные элементы определялись параметрами = 0.1 = 0.1. КИН определялись для сечения у = 0, соответствующих) середине трещины. Результаты сравнения приведены в табл. 2.

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

а б в

Рис. 4. Раскрытие берегов основной и дополнительной трещин в сечении у = 0

В настоящей работе были проведены расчеты в трех случаях расположения двух трещин, соответствующего рис. 1. Во всех трех случаях основная трещина имела радиус, равный п = 1, а дополнительная (малая) трещина — радиус Г2 = 0.5. Начальное расстояние между плоскостями трещин было одинаковым: Н = 0.35. Напряжение, действующее на бесконечности по нормали к плоскости трещин, имело величину — = 0.2. На рис. 4 приведены расчетные положения верхних) и нижних) берегов трещин для трех случаев взаимного расположения (сплошной линией показан верхний берег трещины, пунктирной нижний). Когда малая трещина находится в тени основной, ее влиянием на большую трещину можно практически пренебречь (рис. 4, а и б). В этих случаях раскрытие малой трещины незначительно. Наибольшее влияние малая трещина оказывает в случае, показанном на рис. 4, в: малая трещина раскрыта. Для количественного сравнения влияния трещин друг на друга определялась величина отношения КИН основной трещины к его теоретическому

значению, которое определялось в том случае, когда дополнительная (малая) трещина отсутствует. Координаты ж = — 1 ж ж = 1 соответствуют

у=0

Рис. 3. Положение верхнего берега трещины радиуса Д = 20 с прямоугольными граничными элементами Нхх. Нх = ку = 1

Таблица!

м1/2 Д, м

20 25

Численное значение 0.0063 0.0074

Теоретическое значение 0.0051 0.0056

Т а б л и ц а 2

АУ/у., м1/2 В точке х = — 1 В точке х = 1

Расчет 0.00053 0.00149

Теория 0.00048 0.00139

Т а б л и ц а 3

Случай а Случай б Случай в

х = -1 х = 1 х = -1 х = 1 х = -1 х = 1

0.997 0.998 0.995 0.932 0.999 0.913

Отношение расчетного значения КИН к его аналитическому значению для одиночной трещины (Ki/KT) представлено в табл. 3.

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

Работа выполнена при поддержке гранта РФФИ № 18-01-00151 А.

СПИСОК ЛИТЕРАТУРЫ

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

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

3. Слепян Л.И. Механика трещин. Л.: Судостроение, 1990.

4. Almansi Е. Sull'integrazione dell'equazione differenzigle // Ann. Mat. 1899. 2, N I. 1-51.

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

6. Кузьмин Ю.Н. Осесимметричная деформация упругого слоя, содержащего соосные круговые щели // Изв. АН СССР. Механ. твердого тела. 1972. № 6. 121-130.

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

8. Kassir М.К., Sih G.C. External crack in elastic solid // Int. J. Fract. Mech. 1968. 4, N 4. 347-356.

Поступила в редакцию 21.02.2018