CC BY
22
2018 Математика и механика № 54
УДК 539.375
Б01 10.17223/19988621/54/8
А.В. Малик, И.М. Лавит
МЕТОД РАСЧЕТА КОЭФФИЦИЕНТА ИНТЕНСИВНОСТИ
НАПРЯЖЕНИЙ ДЛЯ НЕПОДВИЖНОЙ ТРЕЩИНЫ НОРМАЛЬНОГО РАЗРЫВА ПРИ ДИНАМИЧЕСКОМ НАГРУЖЕНИИ1
Предлагаемый метод представляет собой модификацию метода прямых. Интегрирование по времени выполняется методом конечных разностей. Краевые задачи, которые получаются на каждом шаге интегрирования по времени, решаются методом конечных элементов. Используются когезионные конечные элементы, обеспечивающие отсутствие сингулярности напряжений в кончике трещины. Результаты расчетов сопоставляются с экспериментальными данными и результатами других исследователей.
Ключевые слова: динамическая механика разрушения, метод прямых, когезионные конечные элементы.
Главным результатом экспериментального исследования динамического разрушения является определение зависимости К() [1], где К1 - коэффициент интенсивности напряжений, t - время. Поэтому сопоставление теории динамического разрушения с экспериментом - это, фактически, сопоставление расчетной и экспериментальной зависимостей Кф.
Методам теоретического определения функции К(() посвящено немало исследований. Аналитические методы ее вычисления [2] применимы только для бесконечных областей. Численные методы представляют собой, как правило, какую-либо модификацию метода конечных элементов или метода граничных элементов. Известно, что высокая точность вычислений достигается при использовании специальных элементов, моделирующих сингулярность напряжений в окрестности кончика трещины [3]. Это усовершенствование сопряжено, однако, с рядом недостатков. Во-первых, сингулярное поле напряжений в окрестности кончика трещины возникает, как правило, через некоторое время после начала процесса динамического нагружения. В течение этого времени сингулярный элемент вынужденно моделирует то, чего на самом деле нет. Это, безусловно, является источником погрешности, которую невозможно оценить заранее.
Второй недостаток обусловлен отличием сингулярных элементов от окружающих их обычных элементов. Для того чтобы сингулярный элемент не вносил погрешность в решение задачи, он должен обеспечивать межэлементную непрерывность перемещений [4]. Такой элемент, моделирующий корневую особенность напряжений в кончике трещины и совместный с четырехугольными изопа-раметрическими элементами второго порядка [4], был сконструирован и успешно применен Барсоумом [5] к решению статических задач. Конечные элементы Бар-соума применены к решению задач динамической механики разрушения в работе [6]. Показано, что их использование позволяет получить достаточно точное реше-
1 Результаты исследований опубликованы при финансовой поддержке ТулГУ в рамках научного проекта № 2017-38ПУБЛ.
ние. Однако подход работы [6] сложен, поэтому в других исследованиях ([3, 7-9] и др.), ориентированных на исследование не только неподвижных, но и растущих трещин, использовались более простые несовместные элементы. Это обусловлено тем, что при моделировании роста трещины с использованием сингулярных элементов приходится перестраивать конечно-элементную сетку. Этот процесс предполагает многократную интерполяцию, что снижает точность расчетов. Погрешность решения, очевидно, тем выше, чем сложнее элементы.
Предложены также методы [10-12], общей чертой которых является отсутствие требования совпадения траектории трещины с границами конечных элементов, а ее вершины - с каким-либо узлом конечно-элементной сетки. Неясно, как можно обосновать сходимость получаемых решений к точным с увеличением числа варьируемых параметров, так как, в отличие от метода конечных элементов, здесь нет непрерывности перемещений и невозможно удовлетворить главным граничным условиям задачи.
Наряду с численными методами, позволяющими непосредственно вычислять коэффициенты интенсивности напряжений, в динамике разрушения используются методы, основанные на представлении о силах сцепления. Первый такой метод был предложен Ксу и Нидлманом [13] и применен к проблеме математического моделирования ветвления трещин. О дальнейшем развитии этого направления можно узнать из обзоров [14-16]. В этих методах удельная высвобожденная энергия представляется как работа сил сцепления. Эту величину можно связать с коэффициентом интенсивности напряжений, если распределение сил сцепления подчинить условиям отсутствия сингулярности поля напряжений в кончике трещины и малости длины зоны сцепления по сравнению с длиной трещины [17, 18]. В упомянутых методах эти условия не выполняются. Исключение составляет метод работы [19]. В нем отсутствие сингулярности удовлетворяется с помощью итерационной процедуры, что является существенным ограничением. Метод был применен [19] лишь к решению квазистатических задач при однопараметриче-ском (простом) нагружении. Решения этим методом задач динамической механики разрушения неизвестны.
Представляет интерес разработка метода, который по точности сопоставим с методами, использующими сингулярные элементы, и обладает простотой и универсальностью методов, основанных на концепции сил сцепления. Такой метод представлен ниже. Он был разработан [20] для решения квазистатических упру-гопластических задач механики разрушения. Его приложение к решению динамических задач не вызвало принципиальных затруднений. Суть метода состоит в использовании специальных когезионных конечных элементов, которые автоматически обеспечивают отсутствие сингулярности напряжений в кончике трещины. Эти элементы таковы, что если не накладывать ограничений на их степени свободы, силы сцепления в пределах элемента будут равны нулю. Поэтому эти элементы можно использовать для расчета напряженного состояния не только в окрестности кончика трещины, но и для всего тела. В этом случае метод приобретает общность, позволяющую решать задачи, рассмотренные Ксу и Нидлманом [13].
Отличие предлагаемого метода от предшествующих заключается также в способе дискретизации задачи. При обычном подходе к решению динамических задач методом конечных элементов перемещения узлов конечно-элементной сетки предполагаются зависящими от времени. В результате после интегрирования по пространственным переменным получается система обыкновенных дифференциальных уравнений относительно узловых неизвестных. В настоящем исследова-
нии использован другой подход, называемый методом прямых [21]. Преобразование континуальной задачи в дискретную осуществляется в два этапа. Вначале исходная начально-краевая задача преобразуется в конечно-разностную задачу по времени. В результате на каждом шаге интегрирования по времени получается краевая задача, решение которой находится методом конечных элементов. Матрица коэффициентов разрешающей системы линейных алгебраических уравнений для всех этих краевых задач одна и та же (задачи отличаются одна от другой за счет правых частей разрешающей системы уравнений), и эта матрица - ленточная. Это свойство удается использовать только в методе прямых.
Постановка задачи
Рассматривается плоская деформация линейно упругого тела. Поперечное сечение тела Б, ограниченное контуром Ь, изображено на рис. 1. Сечение ослаблено прямолинейной трещиной длиной а. К контуру сечения приложена распределенная нагрузка р(*). Сечение и нагрузка симметричны относительно оси абсцисс. Таким образом, исследуемая трещина относится к типу трещин нормального разрыва. Так как ось абсцисс является осью симметрии задачи, ниже рассматривается только верхняя часть поперечного сечения.
Рис. 1. Поперечное сечение тела Fig. 1. Cross-section of the body
Математическая модель основана на соотношениях теории упругости [22]:
Чш = (дкит + дтик )/2; °кш = Щг §кж + 2^кш
дк = д/дхк , X = vE|[(1 + V) (1 - 2v)], ц = Е/[2 (1 + v)], г, к, т = 1,2, (1)
где пк (хт, *) - поле перемещений; &кт - тензор деформаций; окт - тензор напряжений; X, ц - параметры Ламе; Е - модуль Юнга; V - коэффициент Пуассона; дкт - символ Кронекера; р - плотность. В области Б справедливо уравнение движения [22]
дт^кт =РдЛ; Ч =д*ик ; дt =д/д* , (2)
где ук - поле скоростей. Решение уравнения (2) должно удовлетворять начальным
и граничным условиям. Начальные условия предполагаются нулевыми. Граничные условия могут быть либо главными (кинематическими), либо естественными (динамическими) [23]. Трещина считается неподвижной. Требуется определить зависимость К1 (/).
Метод прямых
Слабое решение задачи получается из принципа возможных перемещений
(Ъпк + ) = | РкЬик<ЛЬ , (3)
5 Ь
где 5 - символ вариации. Уравнение (3) эквивалентно уравнению движения (2) и естественным граничным условиям. Решение уравнения (3) разыскивается в классе функций, удовлетворяющих главным граничным условиям. При этом вариации перемещений на участках граничного контура, где заданы главные граничные условия, равны нулю.
Первый шаг решения задачи методом прямых - преобразование уравнения (3) к конечно-разностному уравнению по времени. Для этого в настоящем исследовании используется неявная схема Кранка - Николсон [4]. Пусть Д/ - величина шага интегрирования по времени, п - номер шага интегрирования. Конечно-разностное представление производной по времени на п-м шаге имеет вид
= (Уп - Уп1 )/Д/ ,
(4)
где уп- , уп - значения у на границах временного интервала. Величины, не содержащие производных по времени, представляются на п-м шаге интегрирования по времени в виде
у = ( уп + уп-1 )/2.
Уравнение (3) с учетом выражения для скорости (2) преобразуется к виду
п-1
Д/
5ик + 2(ш + <)) = 2\(р"к + рГ1 )<й,
1 (к + Ук"—1 )
(5)
(6)
2 V - - / Д/
Величины с индексами п — 1 известны из решения для предыдущего шага. Система (6) распадается на два уравнения:
4р
1_(Д/)
2 и"к 5ик +<ш 5Чш
= \( + р"п—1 )5ик<1-
4р
Д/
( и п—1 л к п—1 + VI 1
Д/
5ик—а"к„!5&кш
:
(7)
И 2 {„.п „,п—1\ ,.п—1
Ук = Д (к — ик ГУк .
(8)
Вначале из уравнения (7) определяются перемещения ^, затем из уравнения (8)
- скорости уп
На каждом шаге интегрирования по времени уравнение (7) решается методом конечных элементов. Легко видеть, что и матрица жесткости, и матрица масс на каждом шаге интегрирования по времени будет одна и та же. Пусть и - матрица-столбец узловых неизвестных, V - матрица-столбец их скоростей, Е - матрица жесткости, М -матрица масс, Р - матрица-столбец нагрузок, обусловленная первым слагаемым в правой части уравнения (7). Система линейных алгебраических
уравнений относительно неизвестных ип , эквивалентная вариационному уравнению (7), записывается в виде
.(А' )
-M + E
Un = Pn + Pn-1 +-
(At)
M (U n-1 + AtV n-1)- EU n-1.
Узловые скорости определяются формулой, следующей из равенства (8):
= А (иn
At
U
n-П
тП-1
(9)
(10)
Если пространственное распределение нагрузки неизменно, то уравнение (9) еще больше упрощается, так как можно записать
р п+Р п-1 =[ / (*п)+/ (*п )] рх, (11)
где / (*) - заданная функция - параметр нагружения, Р1 - матрица-столбец нагрузок при / = 1. Ее, так же как и матрицы жесткости и масс, необходимо вычислять только один раз.
Конечные элементы
Сетка конечных элементов состоит из элементов двух типов. Первый из них -это обычные изопараметрические элементы [4]. Второй тип - это специальные ко-гезионные элементы [20]. Они позволяют учитывать специфику, вносимую трещиной в краевые задачи теории упругости. К линии трещины примыкает слой ко-гезионных конечных элементов (рис. 2). Остальная часть сечения покрывается сеткой изопараметрических конечных элементов.
Трещина
Ось симметрии
Рис. 2. Сетка конечных элементов (принципиальная схема). 1 - изопараметрические элементы, 2 - когезионные элементы Fig. 2. Finite element mesh (basic scheme). 1, isoparametric elements and 2, cohesive elements
x
2
В данном исследовании использовались изопараметрические конечные элементы первого порядка. Каждый элемент в локальных координатах представляет собой квадрат (рис. 3).
12
-1 0 -1
11
22
21
Рис. 3. Конечный элемент в локальных координатах. Узлы имеют двойную нумерацию Fig. 3. Finite element in the local coordinates. The nodes have double indexing
П
1
Глобальные координаты точек элемента и перемещения в этих точках определяются формулами [4]:
хк = 11 (5)^ (л)хк,и; ик = 11 (5)^ (п)ики;
1,3 = 1,2; 5,п е [-1,1], (12)
где Хк ц - глобальные координаты узла с номером 13, ик ц - перемещения узла с номером 13, Ь1 (г) - интерполяционные полиномы Лагранжа, которые в данном случае имеют вид
А 00 = 0.5 (1 - г); 12 (г) = 0.5 (1 + г). (13)
Глобальные узловые координаты задаются. Узловые перемещения Пк и - это основные неизвестные задачи.
Предполагается, что в окрестности кончика трещины существует зона сцепления длиной Д, в пределах которой противоположные кромки трещины притягиваются друг к другу [17]. Силы этого притяжения называются силами сцепления. На рис. 4 изображена верхняя кромка трещины в деформированном состоянии. Картина для нижней кромки симметрична.
Излагаемая теория строится на основе постулатов Баренблатта. Согласно первому из них, распределение сил сцепления обеспечивает отсутствие сингулярности поля напряжений в окрестности кончика трещины. Это эквивалентно требованию плавного смыкания кромок трещины в ее кончике (рис. 4), которое записывается в виде
х1 = а, х2 = 0: 51м2 = 0. (14)
Второй постулат утверждает, что длина зоны сцепления намного меньше длины трещины. Третий постулат предполагает, что при продвижении кончика трещины напряженное состояние его ближайшей окрестности не изменяется. Из этих постулатов следует, что и при квазистатическом [17], и при динамическом [18]
нагружении удельная высвобожденная энергия при продвижении кончика трещины и ее связь с коэффициентом интенсивности напряжений определяется формулой [2]:
С 1 -V2
G = 2 J q2diu2dxi = —— К
E
(15)
При выполнении постулатов Баренблатта величина О не зависит от закона изменения сил сцепления вдоль кромки трещины q2 (х1) [17, 18]. Предполагается, что в начальный момент роста трещины скорость движения ее кончика близка к нулю.
Рис. 4. Верхняя кромка трещины при деформации Fig. 4. Top edge of the crack under deformation
Когезионные конечные элементы обеспечивают плавное смыкание кромок трещины в ее кончике и, тем самым, отсутствие сингулярности полей напряжений и деформаций. Когезионные элементы не являются изопараметрическими, и в глобальных координатах представляют собой прямоугольники. При этом формулы связи локальных и глобальных координат упрощаются:
x = ai + М; x2 = a2 + М; ak = 0.5(xk,22 + xkд!);
bk = °.5(Xk,22 - Xk,11 ); X1J1 = X1,12 ; X2,1J = X2
(16)
Перемещения точек когезионного конечного элемента определяются формулами [20]:
«1 = ¿I 00ьи (П)и1,и; «2 = НI ф¿1 (п)и2,л + ЬI фЬ (п)и2Л2 + 0.5^1+2 (п)из,л, 1 ) где Дх - размер элемента по оси абсцисс, из п - значения производной д1и2 в
11-м узле (эти величины определяются только для узлов, лежащих на оси абсцисс - прямой, на которой расположена трещина), Нт (г) - интерполяционные полиномы Эрмита, определяемые формулами [4]:
Н (г) = 0.25(2 - 3г + г3); Н2 (г) = 0.25(2 + 3г - г3);
H3 (z) = 0.25(1 - z - z2 + z3); H4 (z) = 0.25(-1 - z + z2 + z3) .
(18)
Кончик трещины обязательно должен совпадать с каким-либо узлом конечно-элементной сетки. В этом узле полагается не только и2 = 0, но и д«2 = 0, что
a
обеспечивает плавность смыкания кромок трещины в ее кончике. Если не накладывать ограничений на узловые перемещения когезионного конечного элемента, то никаких сил сцепления на его контуре не возникнет. Поэтому область действия сил сцепления (зона сцепления) локализована в пределах конечного элемента, прилегающего к кончику трещины. Таким образом, чем мельче сетка конечных элементов, тем точнее удовлетворяется требование теории Баренблатта о малости длины зоны сцепления по сравнению с длиной трещины.
Далее применяется обычная конечно-элементная процедура, позволяющая свести решение поставленной задачи к решению системы линейных алгебраических уравнений [4]. Главные (кинематические) граничные условия (в рассматриваемом классе задач они могут быть только нулевыми) учитываются следующим образом [4]. Пусть ] - номер узлового перемещения (в глобальной нумерации), которое должно быть равно нулю. Обнуляются ]-я строка и ]-й столбец матрицы коэффициентов разрешающей системы уравнений (последнее - для сохранения симметрии матрицы) за исключением диагонального элемента; обнуляется также }-й элемент столбца свободных членов.
Для вычисления удельной высвобожденной энергии при продвижении трещины используется более простой, чем в работе [20], и, как показывают сравнительные расчеты, более точный подход. Упомянутые граничные условия в кончике трещины представляют собой, с механической точки зрения, связи. Реакции этих связей, QB и MB (рис. 5) легко подсчитать, подставляя значения узловых перемещений в соответствующие непреобразованные уравнения разрешающей системы. При продвижении трещины на длину конечного элемента Дх1 эти реакции связей уменьшаются (по модулю) до нуля, а соответствующие им перемещения ^ и угол поворота д1u2 приобретают конечные значения.
Рис. 5. Деформирование окрестности кончика трещины. Кончик трещины расположен в узле B. Axj - расстояние между соседними узлами, расположенными на линии трещины Fig. 5. Deformation of the crack tip. The crack tip is located in node B. Axj is the distance between adjacent nodes arranged along the crack line
Согласно третьему постулату Баренблатта, при росте трещины конфигурация концевой области не изменяется. Это значит, что перемещения, которые получит узел B при продвижении трещины на величину Atj , будут равны перемещениям
узла A при совпадении кончика трещины с узлом B. При этом величина удельной высвобожденной энергии определяется очевидной формулой:
G = --Д-(0^2Л +MBд^2Л ). (19)
Формула (19) учитывает наличие у трещины двух кромок. Коэффициент интенсивности напряжений вычисляется по формуле, следующей из выражения (15):
к, ^^ Ж
Если из формулы (19) следует, что G < 0, то это значит, что корневая асимптотика еще не установилась и следует считать К, = 0 . Знак минус в формуле (20) соответствует случаю, когда и2Л < 0 . При этом трещина не обязательно закрывается и теория становится неприменимой. Разрез, не имеющий ширины, как модель трещины - это удобная абстракция. В реальных трещинах расстояние между кромками всегда конечно. Отрицательное значение и2Л приводит к уменьшению этого расстояния, но не обязательно к смыканию кромок.
Необходимость введения целого ряда когезионных элементов обусловлена требованием межэлементной непрерывности перемещений, как между когезион-ными элементами, так и между когезионными и обычными изопараметрическими элементами [20].
Численные примеры
Сходимость численных решений исследовалась на различных конечно-элементных сетках при различном количестве шагов по времени N. Относительная погрешность приведенных ниже результатов не превышает 3 %.
Пусть рассматриваемое тело - полоса с центральной трещиной, находящаяся в состоянии плоской деформации. Поперечное сечение полосы изображено на рис. 6. В начальный момент времени к двум противоположным сторонам прямоугольника, параллельным трещине, мгновенно прикладывается равномерно распределенная нагрузка величиной д, далее остающаяся неизменной. Задача имеет две оси симметрии, проходящие через середины противоположных сторон прямоугольника. Поэтому расчетная схема - это четверть сечения, вырезанная осями симметрии (рис. 7).
Граничные условия записываются следующим образом. На участке контура ЛВ: и1 = 0, р2 = 0 - условия симметрии; на участке ВС: р1 = 0, р2 = д; на участке СБ: р1 = 0, р2 = 0; на участке БЕ: и2 = 0, р1 = 0 - условия симметрии; на участке ЕЛ (кромка трещины): р1 = 0, р2 = 0.
Эта задача решалась в работах [6, 24, 25] при следующих исходных данных: V = 0.286 , Н/Ш = 0.385 , а/Ш = 0.23 . Рассматриваемый интервал времени равен tm¡x = 20 мкс. Результаты расчета разработанным методом при N = 220, п1 = 100,
п2 = 40 (квадратная сетка), сопоставляемые с результатами других исследователей [6, 24, 25], приведены на рис. 8. Можно сделать вывод об удовлетворительном согласовании решений всеми четырьмя методами.
2H
'Mttt j
H \
чЧ A / \ W - a 2W / \
4 /
1шн J
q
q
Рис. 6. Прямоугольная область Рис. 7. Расчетная схема
с центральной трещиной Fig. 7. Computational domain
Fig. 6. Rectangular domain with a central crack
sq 1
1 1 1 ■
Ii • 1 2
▲ 3 .
1 ■ 1 1 4
10 t, мкс
15
20
Рис. 8. Результаты решения задачи для прямоугольной области с центральной трещиной. 1 - разработанный метод, 2 - данные работы [24], 3 - данные работы [25], 4 - данные работы [6]
Fig. 8. Solution of the problem for a rectangular domain with a central crack. 1, developed method, 2, data from [24], 3, data from [25], and 4, data from [6]
3
2
0
1
0
5
Рассмотрим теперь задачу о полосе с краевой трещиной. Поперечное сечение полосы изображено на рис. 9. Задача имеет ось симметрии, поэтому расчетная схема - это половина сечения (рис. 7).
2H
H
ш
W
q
q
a
Рис. 9. Прямоугольная область с краевой трещиной Fig. 9. Rectangular domain with an edge crack
Граничные условия записываются следующим образом. На участках контура AB, BC и CD: p1 = 0, p2 = 0; на участке DE: u2 = 0, p1 = 0 - условия симметрии; на участке EA (кромка трещины): pi = 0, p2 = q. Такая расчетная схема моделирует эксперименты Рави-Чандара и Кнаусса [26]. Экспериментальные образцы изготавливались из оптически прозрачного полимера Homalite-100. Его механические характеристики [i]: E = 4550 МПа, v = 0.3i, р = 1230 кг/м3. Размеры образцов:
W = 500 мм, H = 150 мм, a¡W = 0.5. Нагрузка q(t) вначале возрастает по линейному закону до значения qM (при t = tM) и далее остается постоянной. Помимо динамической нагрузки образцы нагружались квазистатически, так как было необходимо несколько раскрыть кромки трещины, чтобы разместить между ними электроды, создающие динамическую нагрузку. В силу линейности задачи коэффициент интенсивности напряжений от статической нагрузки просто прибавлялся к коэффициенту интенсивности напряжений, обусловленному динамической нагрузкой. Расчеты проводились при n1 = 100, n2 = 30. Шаг по времени составлял
At = 10-3 wV^E . Рассмотрены четыре расчетных случая, отличающиеся величиной qM и величиной временного интервала tmax, которая соответствовала экспериментально определенному моменту начала роста трещины (таблица). Величина tM = 25 мкс для всех расчетных случаев.
Максимальные значения нагрузки и временные интервалы [27]
Расчетные случаи qM , МПа tmax > мкс
Рис. 10, a 0.63 102
Рис. 10, b 1.10 56
Рис. 10, с 5.55 18
Рис. 10, d 10.4 17
Рис. 10. Результаты решения задачи и экспериментальные данные для прямоугольной области с краевой трещиной. 1 - разработанный метод, 2 - эксперименты Рави-Чандара и Кнаусса [26, 27]
Fig. 10. Problem solution and experimental data for a rectangular domain with an edge crack. 1, developed method and 2, experiments by Ravi-Chandar and Knauss [26, 27]
Сопоставление результатов расчетов с экспериментальными данными приведено на рис. 10. Можно сделать вывод об удовлетворительном согласовании теории с экспериментом.
Заключение
Предложен новый метод решения задач динамической механики разрушения. Метод позволяет вычислять зависимость коэффициента интенсивности напряжений от времени для плоскодеформированного тела с неподвижной трещиной нормального разрыва. Суть метода состоит в использовании специальных когезион-ных элементов [20]. Эти элементы позволяют получить решение задачи механики разрушения для модели когезионной трещины Баренблатта. По силам реакции связей и узловым параметрам находится удельная высвобожденная энергия и по ней, вследствие эквивалентности теорий Гриффитса и Баренблатта [18], определяется коэффициент интенсивности напряжений. Когезионные элементы удовлетворяют требованиям межэлементной непрерывности перемещений и полноты. Другая важная особенность предложенного метода - это получение решения динамической задачи методом прямых.
Возможности метода продемонстрированы на решении двух задач для прямоугольных областей с трещинами. Решения сопоставлялись с решениями других авторов и с экспериментальными данными. Можно отметить их удовлетворительное согласование.
К достоинствам метода можно отнести также его простоту и отсутствие необходимости перестраивать конечно-элементную сетку при росте трещины. Это дает принципиальную возможность использовать метод для математического моделирования роста трещин.
ЛИТЕРАТУРА
1. Ravi-Chandar K. Dynamic fracture. Oxford: Elsevier, 2004.
2. FreundL.B. Dynamic fracture mechanics. New York: Cambridge University Press, 1998.
3. Нисиока Т., Атлури C. Вычислительные методы в динамике разрушения // Вычислительные методы в механике разрушения. М.: Мир, 1990. С. 267-321.
4. Zienkiewicz O.C., Taylor R.L. The finite element method. Oxford: Butterworth-Heinemann, 2000.
5. Barsoum R.S. Application of quadratic isoparametric finite elements in linear fracture // Int. J. Fract. 1974. V. 10. P. 603-605.
6. Murti V., Valliappan S. The use of quarter point element in dynamic crack analysis // Eng. Fract. Mech. 1986. V. 23. P. 585-614.
7. Liu P., Bui T.Q., Zhang C., Yu T.T., Liu G.R., Golub M.V. The singular edge-based smoothed finite element method for stationary dynamic crack problems in 2D elastic solids // Comput. Methods Appl. Mech. Engrg. 2012. V. 233-236. P. 68-80.
8. Wu L., Liu P., Shi C., Zhang Z., Bui T.Q., Jiao D. Edge-based smoothed extended finite element method for dynamic fracture analysis // Appl. Math. Modelling. 2016. V. 40. P. 8564-8579.
9. Yao W., Cai Z., Hu X. A new symplectic analytical singular element for crack problems under dynamic loading condition // Eng. Fract. Mech. 2018. V. 188. P. 431-447.
10. Menouillard T., Belytschko T. Dynamic fracture with meshfree enriched XFEM // Acta Mech. 2010. V. 213. P. 53-69.
11. Kumar S., Singh I.V., Mishra B.K., Singh A. New enrichments in XFEM to model dynamic crack response of 2-D elastic solids // Int. J. Impact Engrg. 2016. V. 87. P. 198-211.
12. Li W., Zheng H., Sun G. The moving least squares based numerical manifold method for vibration and impact analysis of cracked bodies // Eng. Fract. Mech. 2018. V. 190. P. 410434.
13. Xu X.P., Needleman A. Numerical simulations of fast crack growth in brittle solids // J. Mech. Phys. Solids. 1994. V. 42. P. 1397-1434.
14. de Borst R., Remmers J.J.C., Needleman A. Mesh-independent discrete numerical representations of cohesive-zone models // Eng. Fract. Mech. 2006. V. 73. P. 160-177.
15. Song J.H., Wang H., Belytschko T. A comparative study on finite element methods for dynamic fracture // Comput. Mech. 2008. V. 42. P. 239-250.
16. Sukumar N., Dolbow J.E., Moes N. Extended finite element method in computational fracture mechanics: a retrospective examination // Int. J. Fract. 2015. V. 196. P. 189-206.
17. Баренблатт Г.И. Математическая теория равновесных трещин, образующихся при хрупком разрушении // Журн. прикл. мех-ки и техн. физики. 1961. № 4. С. 3-56.
18. Willis J.R. A comparison of the fracture criteria of Griffith and Barenblatt // J. Mech. Phys. Solids. 1967. V. 15. P. 151-162.
19. Moes N., Belytschko T. Extended finite element method for cohesive crack growth // Eng. Fract. Mech. 2002. V. 69. P. 813-833.
20. Лавит И.М. Об устойчивом росте трещины в упругопластическом материале // Проблемы прочности. 1988. № 7. С. 18-23.
21. Ладыженская О.А. Краевые задачи математической физики. М.: Наука, 1973. 408 с.
22. Тимошенко С.П., Гудьер Дж. Теория упругости. М.: Наука, 1975. 576 с.
23. Михлин С.Г. Вариационные методы в математической физике. М.: Наука, 1970. 512 с.
24. Nishioka T., Atluri S.N. Numerical modeling of dynamic crack propagation in finite bodies, by moving singular elements. Part 2: Results // J. Appl. Mech. 1980. V. 47. P. 577-582.
25. Fedelinski P., Aliabadi M.H., Rooke D.P. The dual boundary element method: J-integral for dynamic stress intensity factor // Int. J. Fract. 1994. V. 65. P. 369-381.
26. Ravi-Chandar K., Knauss W.G. Dynamic crack-tip stresses under stress wave loading - a comparison of theory and experiment // Int. J. Fract. 1982. V. 20. P. 209-222.
27. Ma C.C., FreundL.B. The extent of the stress intensity factor field during crack growth under dynamic loading conditions // J. Appl. Mech. 1986. V. 53. P. 303-310.
Статья поступила 20.03.2018 г.
Malik A.V., Lavit I.M. (2018) ON THE COMPUTATION METHOD FOR THE STRESS INTENSITY FACTOR OF A STATIONARY CRACK IN MODE I UNDER DYNAMIC LOADING. Vestnik Tomskogo gosudarstvennogo universiteta. Matematika i mekhanika [Tomsk State University Journal of Mathematics and Mechanics]. 54. pp. 88-102
DOI 10.17223/19988621/54/8
Keywords: dynamic fracture, method of lines, cohesive finite elements.
A method for calculating the time dependence of the stress intensity factor of a dynamic loaded body with a stationary crack in mode I is proposed. The method represents a modified method of lines developed for solving the problems of dynamic fracture mechanics. Time integration is performed using the finite difference method. The Crank-Nicolson implicit scheme is applied. Boundary problems obtained at each step of time integration are solved by the finite element method. For each time instant, the stress intensity factor is determined by the calculated value of specific energy release. For this purpose, the special cohesive finite elements ensured by the implementation of Barenblatt's postulates are used. Addition of the degrees of freedom for the mesh points arranged along the crack line enables to provide a smooth joining of the edges at the crack tip. This is equivalent to the absence of stress singularity in the tip of crack. The results of calculations are compared with those obtained by other researchers and with experimental data. A satisfactory agreement ensures the efficiency of the method applied. The latter can also be used to solve the problems of growing and branching cracks. Moreover, it admits taking into account plastic deformation.
MALIK Alexander Vasil'evich (Tula State University, Tula, Russian Federation). E-mail: [email protected]
LAVIT Igor Mikhaylovich (Doctor of Physics and Mathematics, Tula State University, Tula, Russian Federation). E-mail: [email protected]
REFERENCES
1. Ravi-Chandar K. (2004) Dynamic Fracture. Oxford: Elsevier.
2. Freund L.B. (1998) Dynamic Fracture Mechanics. New York: Cambridge University Press.
3. Nishioka T., Atluri S.N. (1986) Computational methods in dynamics fracture. Computational Methods in the Mechanics of Fracture. Amsterdam: Elsevier. pp. 335-383.
4. Zienkiewicz O.C., Taylor R.L. (2000) The finite element method. Oxford: ButterworthHeinemann.
5. Barsoum R.S. (1974) Application of quadratic isoparametric finite elements in linear fracture. Int. J. Fract. 10. pp. 603-605. DOI: 10.1007/BF00155266.
6. Murti V., Valliappan S. (1986) The use of quarter point element in dynamic crack analysis. Eng. Fract. Mech. 23. pp. 585-614. DOI: 10.1016/0013-7944(86)90164-5.
7. Liu P., Bui T.Q., Zhang C., Yu T.T., Liu G.R., Golub M.V. (2012) The singular edge-based smoothed finite element method for stationary dynamic crack problems in 2D elastic solids. Comput. Methods Appl. Mech. Engrg. 233-236. pp. 68-80. DOI: 10.1016/j.cma.2012.04.008.
102
А.В. Manm, M.M. flaBUT
8. Wu L., Liu P., Shi C., Zhang Z., Bui T.Q., Jiao D. (2016) Edge-based smoothed extended finite element method for dynamic fracture analysis. Appl. Math. Modelling. 40. pp. 85648579. DOI: 10.1016/j.apm.2016.05.027.
9. Yao W., Cai Z., Hu X. (2018) A new symplectic analytical singular element for crack problems under dynamic loading condition. Eng. Fract. Mech. 188. pp. 431-447. DOI: 10.1016/j.engfracmech.2017.09.016.
10. Menouillard T., Belytschko T. (2010) Dynamic fracture with meshfree enriched XFEM. Acta Mech. 213. pp. 53-69. DOI: 10.1007/s00707-009-0275-z.
11. Kumar S., Singh I.V., Mishra B.K., Singh A. (2016) New enrichments in XFEM to model dynamic crack response of 2-D elastic solids. Int. J. Impact Eng. 87. pp. 198-211. DOI: 10.1016/j.ijimpeng.2015.03.005.
12. Li W., Zheng H., Sun G. (2018) The moving least squares based numerical manifold method for vibration and impact analysis of cracked bodies. Eng. Fract. Mech. 190. pp. 410-434. DOI: 10.1016/j.engfracmech.2017.12.025.
13. Xu X.P., Needleman A. (1994) Numerical simulations of fast crack growth in brittle solids. J. Mech. Phys. Solids. 42. pp. 1397-1434. DOI: 10.1016/0022-5096(94)90003-5.
14. De Borst R., Remmers J.J.C., Needleman A. (2006) Mesh-independent discrete numerical representations of cohesive-zone models. Eng. Fract. Mech. 73. pp. 160-177. DOI: 10.1016/ j.engfracmech.2005.05.007.
15. Song J.H., Wang H., Belytschko T. (2008) A comparative study on finite element methods for dynamic fracture. Comput. Mech. 42. pp. 239-250. DOI: 10.1007/s00466-007-0210-x.
16. Sukumar N., Dolbow J.E., Moës N. (2015) Extended finite element method in computational fracture mechanics: a retrospective examination. Int. J. Fract. 196. pp. 189-206. DOI: 10.1007/s1070.
17. Barenblatt G.I. (1962) The mathematical theory of equilibrium cracks in brittle fracture. Adv. Appl. Mech. 7. pp. 55-129. DOI: 10.1016/S0065-2156(08)70121-2.
18. Willis J.R. (1967) A comparison of the fracture criteria of Griffith and Barenblatt. J. Mech. Phys. Solids. 15. pp. 151-162. DOI: 10.1016/0022-5096(67)90029-4.
19. Moës N., Belytschko T. (2002) Extended finite element method for cohesive crack growth. Eng. Fract. Mech. 69. pp. 813-833. DOI: 10.1016/S0013-7944(01)00128-X.
20. Lavit I.M. (1988) Stable crack growth in an elastoplastic material. Strength of Materials. 20. pp. 854-860. DOI: 10.1007/BF01528695.
21. Ladyzhenskaya O.A. (1985) The Boundary Value Problems of Mathematical Physics. Berlin: Springer-Verlag.
22. Timoshenko S.P., Goodier J.N. (1970) Theory of Elasticity. New York: McGraw-Hill.
23. Mikhlin S.G. (1964) Variational Methods in Mathematical Physics. Oxford: Pergamon Press.
24. Nishioka T., Atluri S.N. (1980) Numerical modeling of dynamic crack propagation in finite bodies, by moving singular elements. Part 2: Results. J. Appl. Mech. 47. pp. 577-582. DOI: 10.1115/1.3153734.
25. Fedelinski P., Aliabadi M.H., Rooke D.P. (1994) The dual boundary element method: J-integral for dynamic stress intensity factor. Int. J. Fract. 65. pp. 369-381. DOI: 10.1007/ BF00012375.
26. Ravi-Chandar K., Knauss W.G. (1982) Dynamic crack-tip stresses under stress wave loading - a comparison of theory and experiment. Int. J. Fract. 20. pp. 209-222. DOI: 10.1007/ BF01140336.
27. Ma C.C., Freund L.B. (1986) The extent of the stress intensity factor field during crack growth under dynamic loading conditions. J. Appl. Mech. 53. pp. 303-310. DOI: 10.1115/ 1.3171756.
