A shock loading on a bar with a central crack

Abstract


Linear thermoelasticity is studied of a plane regular truss formed by four families of The paper is concerned with the problem of calculating the time dependence of the stress intensity factor for a plane-strain bar with a stationary central crack caused by the opening mode. A uniformly distributed load has been immediately imposed on the basis of the bar and remained unchanged later. The model of the crack with cohesive forces distributed by Barenblatt’s postulates is used. In this case the stress intensity factor is a result of the calculation of the released energy that is determined by the cohesive forces. The solution is found with a new numerical method that is an adaptation of the method of lines to dynamic fracture mechanics problems. The Crank-Nicolson implicit finite-difference scheme is used for time integration. Boundary problems arising at each step of time integration are solved by the finite element method. The special cohesive finite elements are used, so that the solution of the problem could satisfy Barenblatt’s postulates. Previously these elements were used to solve quasi-static nonlinear fracture mechanics problems. By introducing the additional degrees of freedom of the nodes lying on the crack line, it becomes possible to ensure a smooth closing of the crack edges at its tip; and it is equivalent to the absence of the singularity stress and strain fields at its tip. The cohesive forces are calculated as the constraints. Their field of action (cohesive zone) is localized within the finite element which is adjacent to the tip of the crack. Thus, the smaller the finite element mesh is, the better it satisfies the requirement of Barenblatt’s theory. This requirement concerns the length of the cohesive zone which is small compared to the length of the crack. The stated problem is called Chen’s problem and had earlier been solved by researchers that used different methods. The proximity of the obtained results makes it possible to consider Chen’s problem as a test; and its solution obtained by the developed method agrees well with the data of other researchers.

Full Text

Введение Формулировка условий старта трещины при динамическом нагружении - сложная проблема, до сих пор окончательно не решенная [1-3]. Хотя предлагаемые теории различны (см., например, работы [4-6]), все они включают в себя как составную часть зависимость коэффициента интенсивности напряжений (КИН) от времени. Ниже везде речь идет только о трещинах нормального разрыва, и поэтому под КИН понимается величина . Определение функции , где - время, представляет собой самостоятельную задачу, аналитическое решение которой возможно только для бесконечных областей [1-3]. Существующие численные методы (см. обзоры [7-10]) обладают теми или иными недостатками, поэтому разработка новых эффективных методов расчета представляет интерес. По-видимому, наиболее перспективным направлением вычислительной динамики разрушения является создание методов, основанных на модели когезионной трещины, кромки которой притягиваются друг к другу силами сцепления [11]. Главное достоинство этих методов - алгоритмическая простота, позволяющая решать ими сложные задачи. Первый такой метод был предложен в работе [12] и применен к исследованию ветвления трещин. Дальнейшее развитие - это, как правило, объединение подхода работы [12] с тем или иным вариантом метода конечных элементов, а также использование различных конституционных соотношений для сил сцепления (см., например, работы [11, 13-16]). Для того чтобы вычислять КИН, используя модель когезионной трещины, необходимо, чтобы распределение сил сцепления подчинялось постулатам Баренблатта [17, 18, 2]. Большинство методов, в том числе и упомянутые выше, этому требованию не удовлетворяет. Одним из исключений является метод, изложенный в работе [19] и предназначенный для решения квазистатических задач. Его обобщение на динамические задачи неизвестно и, по-видимому, невозможно, так как в нем существенным этапом является вычисление -интеграла, а последний перестает быть инвариантным при необходимости учитывать инерционные силы. Есть и другой метод, который также был применен ранее для решения квазистатических задач [20]. Этот метод допускает обобщение на динамические задачи. Оно представлено ниже. В нем, по сравнению с работой [20], учтены силы инерции и зависимость всех полевых величин от времени. Некоторые аспекты получившегося метода решения динамических задач приведены в статье [21]. Он представляет собой модификацию метода прямых (метода Роте) [22-24]. С использованием конечно-разностной дискретизации по времени решение задачи сводится к совокупности последовательно решаемых краевых задач теории упругости. Выбор конечно-разностной схемы интегрирования по времени неоднозначен: возможно применение любой надежной схемы, как явной, так и неявной. В настоящей работе использована широко применяемая, в том числе и совместно с методом конечных элементов, неявная схема Кранка-Николсон [25]. Для решения краевых задач используется метод конечных элементов [25]. Проблема согласования распределения сил сцепления с постулатами Баренблатта решается с помощью включения в конечно-элементную сетку специальных когезионных элементов [20]. Оценка эффективности численного метода основывается на решении тестовых задач. Обычно эти задачи имеют аналитическое решение. Но в механике разрушения таких задач для конечных областей нет, и поэтому тестовые задачи приходится выбирать из задач, решаемых численно. И если подходящих статических задач можно найти немало [26], то для случая внезапно прикладываемой нагрузки - только одну. Это задача о полосе с неподвижной центральной трещиной; к основаниям полосы мгновенно прикладывается равномерно распределенная нагрузка, остающаяся далее неизменной. Впервые эту задачу решил Чен [27]; она была предметом ряда последующих исследований ([1, 28-33] и др.). В работе [30] эта задача названа задачей Чена, и это название закрепилось в более поздних публикациях. Считать задачу Чена тестовой позволяет то, что авторы упомянутых исследований получили близкие результаты, причем они применяли различные методы. Так, в работах [27, 30] использовался метод конечных разностей для интегрирования и по времени, и по пространственным координатам с применением явной трехслойной схемы. В работах [1, 28] для интегрирования по пространству использовался метод конечных элементов, причем в работе [1] для моделирования особенности поля напряжений в кончике трещины применялись специальные сингулярные конечные элементы. В результате получалась система обыкновенных дифференциальных уравнений, которая интегрировалась численно по явной трехслойной схеме Ньюмарка. В работах [29, 31-33] применялся метод граничных элементов. В работе [29] для моделирования особенности поля напряжений использовались сингулярные граничные элементы; были найдены зависящие от времени функции Грина, благодаря чему интегрирование по времени свелось к вычислению определенных интегралов. В работах [31, 32] применялось преобразование Лапласа по времени, а в работе [33] - преобразование Фурье. 1. Постановка задачи Рассматриваются малые деформации изотропного однородного линейно-упругого тела. Тело содержит трещину, длина которой неизменна; форма тела и распределение нагрузок таковы, что эта трещина является трещиной нормального отрыва, то есть где , , - коэффициенты интенсивности напряжений [1]. Пусть рассматриваемое тело - полоса, находящаяся в состоянии плоской деформации. Поперечное сечение полосы изображено на рис. 1. В начальный момент времени к двум противоположным сторонам прямоугольника, параллельным трещине, мгновенно прикладывается равномерно распределенная нагрузка величиной , далее остающаяся неизменной. Ставится задача определить зависимость . Рис. 1. Поперечное сечение полосы - прямоугольник с симметрично расположенной прямолинейной центральной трещиной Fig. 1. The cross-section of the strip. This is the rectangle with a symmetrical rectilinear central crack Задача имеет две оси симметрии, проходящие через середины противоположных сторон прямоугольника. Поэтому расчетная схема - это четверть сечения, вырезанная осями симметрии (рис. 2). Рис. 2. Расчетная схема задачи Чена Fig. 2. The design model of Chen’s problem Граничные условия задачи определяются условиями нагружения и условиями симметрии. Пусть - вектор перемещений в декартовой системе координат (см. рис. 2); - внешняя распределенная нагрузка; . Граничные условия записываются следующим образом. На участке контура : , - условия симметрии; на участке : , ; на участке : , ; на участке : , - условия симметрии; на участке (кромка трещины): , . Деформирование полосы определяется следующими соотношениями [1, 2]: (1) где - тензор деформаций; - символ Кронекера; - средняя деформация, - объемный модуль и модуль сдвига; - тензор напряжений; - плотность материала. В начальный момент времени () полоса покоится. 2. Метод решения задачи Численное решение задачи получается на основе принципа возможных перемещений (2) где - площадь области (см. рис. 1); - ее граничный контур; - символ вариации; - поле скоростей. В соответствии с методом прямых [22-24] перейдем в уравнении (2) к конечно-разностному представлению, используя неявную схему Кранка-Николсон [25]. Пусть - величина шага интегрирования по времени, - номер шага. Конечно-разностное представление производной по времени на п-м шаге имеет вид где - значения на границах -го временного интеграла (шага). Величины, не содержащие производных по времени, представляются на п-м шаге интегрирования по времени в виде С учетом этих формул конечно-разностное представление уравнения (2) записывается в виде (3) Величины с индексами известны из решения для предыдущего шага. Из системы (3) определяются величины с индексами . Выразим из второго уравнения системы (3) величину и подставим в первое уравнение. Получим (4) (5) Таким образом, система (3) распадается на два последовательно решаемых уравнения. Вначале из уравнения (5) определяются перемещения , затем из уравнения (4) - скорости . Уравнение (5) представляет собой вариационное уравнение по пространственным переменным, которое необходимо решать на каждом шаге интегрирования по времени. Его решение находится методом конечных элементов [25]. Специфика, вносимая трещиной в краевую задачу теорию упругости, учитывается с помощью специальных когезионных конечных элементов - так, как это описано в работе [20]. Рис. 3. Конечно-элементная сетка (принципиальная схема). Когезионные элементы заштрихованы Fig. 3. Finite-element mesh (schematic). Cohesive elements are shaded Когезионные элементы составляют горизонтальный ряд, прилегающий к кромке трещины (рис. 3). В локальных координатах [25] и обычные, и когезионные элементы представляют собой одинаковые квадраты. Глобальные координаты точек элемента определяются как где - заданные глобальные декартовы координаты узлов (верхние индексы обозначают номер узла в локальной нумерации); - локальные координаты точек элемента; - интерполяционные полиномы Лагранжа. Перемещения точек обычного конечного элемента задаются аналогичной формулой [25] (6) где - узловые перемещения. Перемещения точек когезионного конечного элемента определяются формулами [20] (7) где - размер элемента по оси абсцисс; - значения производной в ij-м узле (эти величины определяются только для узлов, лежащих на оси абсцисс - прямой, на которой расположена трещина); - интерполяционные полиномы Эрмита, Представление перемещений в виде (6), (7) обеспечивает межэлементную непрерывность поля перемещений во всей области Необходимость введения дополнительных узловых степеней свободы обусловлена требованием плавного смыкания кромок трещины в ее кончике [17]. Конечно-элементная сетка строится так, чтобы кончик трещины совпадал с каким-либо узлом. В этом узле полагается не только , но и , что и обеспечивает упомянутую плавность смыкания кромок трещины. Эти условия представляют собой с механической точки зрения связи, уменьшающие число степеней свободы узла, совпадающего с кончиком трещины. Реакции этих связей, и (рис. 4), определяются обычным образом после решения системы линейных алгебраических уравнений относительно узловых неизвестных. При продвижении трещины на длину конечного элемента эти реакции связей уменьшаются (по модулю) до нуля, а соответствующие им перемещения и угол поворота приобретают конечные значения. Рис. 4. Реакции связей и узловые перемещения в узлах конечного элемента, прилегающего к кончику трещины (узел ) Fig. 4. Reactions of ties and nodal displacements at the nodes of the finite element which is adjacent to the crack tip (node ) Обычно принимается допущение о неизменяемости конфигурации концевой области при продвижении трещины [17]. Это значит, что перемещения, которые получит узел при продвижении трещины на величину , будут равны перемещениям узла при совпадении кончика трещины с узлом . При этом величина высвобожденной энергии, отнесенная к приращению длины трещины , определяется как (8) Здесь учтено, что трещина имеет две кромки. Можно показать [17, 18], что на расстояниях, значительно превышающих длину когезионной зоны, поле напряжений определяется асимптотическими формулами [1-3], получаемыми при решении задач при отсутствии сил сцепления. При этом связь между КИН и величиной дается формулой [1, 2] , (9) где - модуль Юнга; - коэффициент Пуассона. Если из формулы (8) получается, что , то это значит, что корневая асимптотика [1-3] еще не установилась и, следовательно, КИН равен нулю. Знак минус в формуле (9) соответствует случаю, когда . Не следует считать, что при этом трещина закрывается и теория становится неприменимой. Разрез, не имеющий ширины, как модель трещины - это удобная абстракция. В реальных трещинах расстояние между кромками всегда конечно. Отрицательное значение приводит к уменьшению этого расстояния, но не обязательно к смыканию кромок. 3. Решение задачи Чена Расчеты проводились при следующих исходных данных [27]: модуль Юнга МПа, коэффициент Пуассона , плотность кг/м3, нагрузка МПа, размеры полосы: мм, мм, длина трещины мм. Рассматриваемый интервал времени мкс. Объемный модуль и модуль сдвига выражаются через и известными формулами Сходимость численного решения исследовалась на различных конечно-элементных сетках при различном количестве шагов по времени N. В таблице приведены результаты вычисления максимальной величины безразмерного КИН для трех вариантов расчета, отличающихся значениями и числом конечных элементов. Во всех вариантах отношение числа конечных элементов по оси ординат к числу конечных элементов по оси абсцисс принималось постоянным, . Максимальные значения для различных вариантов расчета The highest values of for different calculation variants 450 50 2,639 900 50 2,652 900 100 2,648 На рис. 5 представлены результаты расчета при , и . При этом относительная погрешность расчета, согласно данным таблицы, не превышает 1 %. На рис. 5 по оси абсцисс отложено безразмерное время [30] где - скорость волны расширения [1-3], Рис. 5. Зависимость КИН от времени: 1 - решение Чена [27]; 2 - результаты работы [30]; 3 - результаты расчета разработанным методом Fig. 5. The dependence of stress intensity factor on time: 1 is Chen’s solution [27]; 2 are the results of [30]; 3 are the calculation results based on the developed method Как следует из рис. 5, результаты расчетов изложенным методом практически совпадают с результатами работ [27, 30]. Резкие качественные изменения различных участков графиков обусловлены распространением, взаимодействием и отражением от границ области волн расширения, сдвига и волн Рэлея. Обсуждение этих эффектов можно найти в работах [1, 27, 30]. Заключение Решение задачи Чена изложенным методом свидетельствует о его приемлемой точности. Другие его достоинства: простота и возможность применения к более сложным динамическим задачам, в частности задачам о распространении трещин, в том числе и с учетом пластического деформирования. Примеры решения таких задач в квазистатической постановке даны в работе [20].

About the authors

A V Malik

Tula State University

I E Ryazantseva

Tula State University

I M Lavit

Tula State University

References

  1. Партон В.З., Борисковский В.Г. Динамика хрупкого разрушения. - М.: Машиностроение, 1988. - 240 с.
  2. Freund L.B. Dynamic fracture mechanics. - New York: Cambridge university press, 1998. - 563 p.
  3. Ravi-Chandar K. Dynamic fracture. - Amsterdam: Elsevier, 2004. - 254 p.
  4. Морозов Н.Ф., Петров Ю.В. Проблемы динамики разрушения твердых тел. - СПб.: Изд-во С.-Петерб. у-та, 1997. - 129 с.
  5. Zhao Y.-P. Suggestion of a new criterion of dynamic fracture // Int. J. Fract. - 1995. - Vol. 71. - P. R77-R78.
  6. Liu C., Knauss W.G., Rosakis A.J. On the modeling of fracture of brittle solids // Int. J. Fract. - 1998. - Vol. 90. - P. 103-118.
  7. Song J.-H., Wang H., Belytschko T. A comparative study on finite element methods for dynamic fracture // Comput. Mech. - 2008. - Vol. 42. - P. 239-250.
  8. Братов В.А. Численные модели динамики разрушения // Вычислительная механика сплошных сред. - 2009. - Т. 2, № 3. - С. 5-16.
  9. Rabczuk T. Computational methods for fracture in brittle and quasi-brittle solids: state-of-the-art review and future perspectives // ISRN Applied Mathematics. - 2013. - Vol. 2013. - Art. ID849231. - 38 p. - URL: http://dx.doi.org/10.1155/2013/849231.
  10. Fineberg J., Bouchbinder E. Recent developments in dynamic fracture: some perspectives // Int. J. Fract. - 2015. - Vol. 196. - P. 33-57.
  11. The cohesive zone model: advantages, limitations and challenges / M. Elices, G.V. Guinea, J. Gomez, J. Planas // Eng. Fract. Mech. - 2002. - Vol. 69. - P. 137-163.
  12. Xu X.-P., Needleman A. Numerical simulations of fast crack growth in brittle solids // J. Mech. Phys. Solids. - 1994. - Vol. 42. - P. 1397-1434.
  13. De Borst R., Remmers J.J.C., Needleman A. Mesh-independent discrete numerical representations of cohesive-zone models // Eng. Fract. Mech. - 2006. - Vol. 73. - P. 160-177.
  14. Bardenhagen S.G., Nairn J.A., Lu H. Simulation of dynamic fracture with the material point method using a mixed J-integral and cohesive law approach // Int. J. Fract. - 2011. - Vol. 170. - P. 49-66.
  15. Agwai A., Guven I., Madenci E. Predicting crack propagation with peridynamics: a comparative study // Int. J. Fract. - 2011. - Vol. 171. - P. 65-78.
  16. Javidrad F., Mashayekhy M. A cohesive zone model for crack growth simulation in AISI 304 steel // J. Solid Mechanics. - 2014. - Vol. 6. - P. 378-388.
  17. Баренблатт Г.И. Математическая теория равновесных трещин, образующихся при хрупком разрушении // Журн. прикл. механики и техн. физики. - 1961. - № 4. - С. 3-56.
  18. Willis J.R. A comparison of the fracture criteria of Griffith and Barenblatt // J. Mech. Phys. Solids. - 1967. - Vol. 15. - P. 151-162.
  19. Moes N., Belytschko T. Extended finite element method for cohesive crack growth // Eng. Fract. Mech. - 2002. - Vol. 69. - P. 813-833.
  20. Лавит И.М. Об устойчивом росте трещины в упругопластическом материале // Проблемы прочности. - 1988. - № 7. - С. 18-23.
  21. Малик А.В., Белая Л.А., Лавит И.М. О динамическом нагружении тела с трещиной в условиях плоской деформации // Фундаментальные и прикладные проблемы техники и технологии. - 2016. - № 1 (315). - С. 3-10.
  22. Михлин С.Г. Прямые методы в математической физике. - М.: ГИТТЛ, 1950. - 428 с.
  23. Лисковец О.А. Метод прямых // Дифференциальные уравнения. - 1965. - Т. 1, № 12. - С. 1662-1678.
  24. Ладыженская О.А. Краевые задачи математической физики. - М.: Наука, 1973. - 408 с.
  25. Зенкевич О. Метод конечных элементов в технике. - М.: Мир, 1975. - 541 с.
  26. Справочник по коэффициентам интенсивности напряжений. Т. 1 / под ред. Ю. Мураками. - М.: Мир, 1990. - 448 с.
  27. Chen Y.M. Numerical computation of dynamic stress intensity factors by a lagrangian finite-difference method (the HEMP code) // Eng. Fract. Mech. - 1975. - Vol. 7. - P. 653-660.
  28. Brickstad B. A FEM analysis of crack arrest experiments // Int. J. Fract. - 1983. - Vol. 21. - P. 177-194.
  29. Israil A.S.M., Dargush G.F. Dynamic fracture mechanics studies by time-domain BEM // Eng. Fract. Mech. - 1991. - Vol. 39. - P. 315-328.
  30. Lin X., Ballmann J. Re-consideration of Chen’s problem by finite difference method // Eng. Fract. Mech. - 1993. - Vol. 44. - P. 735-739.
  31. Wen P.H., Aliabadi M.H., Rooke D.P. Application of the weight function method to two-dimensional elastodynamics fracture mechanics // Int. J. Fract. - 1996. - Vol. 76. - P. 193-206.
  32. Wen P.H., Aliabadi M.H., Rooke D.P. A contour integral method for dynamic stress intensity factors // Theor. and Appl. Fract. Mechanics. - 1997. - Vol. 27. - P. 29-41.
  33. Phan A.-V. A non-singular boundary integral formula for frequency domain analysis of the dynamic T-stress // Int. J. Fract. - 2012. - Vol. 173. - P. 37-48.

Statistics

Views

Abstract - 289

PDF (Russian) - 96

Cited-By


PlumX


Copyright (c) 2017 Malik A.V., Ryazantseva I.E., Lavit I.M.

Creative Commons License
This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.

This website uses cookies

You consent to our cookies if you continue to use our website.

About Cookies