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

Научная статья DOI: 10.18287/2541-7525-2023-29-1-47-63

УДК 512.531; 519.7 Дата: поступления статьи: 03.03.2023

после рецензирования: 14.04.2023 принятия статьи: 30.05.2023

Д.В. Чаплий

Самарский национальный исследовательский университет имени академика С.П. Королева, г. Самара, Российская Федерация E-mail: [email protected]. ORCID: https://orcid.org/0000-0001-9510-3659

Л.В. Степанова

Самарский национальный исследовательский университет имени академика С.П. Королева, г. Самара, Российская Федерация E-mail: [email protected] ORCID: https://orcid.org/0000-0002-6693-3132

О.Н. Белова

Самарский национальный исследовательский университет имени академика С.П. Королева, г. Самара, Российская Федерация E-mail: [email protected] ORCID: https://orcid.org/0000-0002-4492-223X

ВОЗДЕЙСТВИЕ АККУМУЛЯЦИИ ПОВРЕЖДЕНИЙ НА АСИМПТОТИЧЕСКОЕ ПОВЕДЕНИЕ НАПРЯЖЕНИЙ В ОКРЕСТНОСТИ ВЕРШИНЫ ТРЕЩИНЫ1

АННОТАЦИЯ

Предметом данного исследования является анализ механических полей, ассоциированных с вершиной трещины, находящейся в условиях ползучести при учете явления накопления повреждений. Задачей работы является проведение конечно-элементного моделирования, с помощью программного комплекса SIMULIA Abaqus, одноосного растяжения пластины с центральной горизонтальной трещиной в условиях ползучести с учетом накопления поврежденности в плоской постановке задачи. Для численного моделирования ползучести используется степенной закон Бейли — Нортона. Степенной закон ползучести с помощью пользовательской процедуры UMAT (User Material) пакета SIMULIA Abaqus был дополнен кинетическим уравнением накопления поврежденности Качанова — Работнова в связанной постановке. В данном исследовании вершина трещины моделировалась в виде математического надреза, а также выреза с конечным радиусом закругления вершин.

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

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

Ключевые слова: пользовательская процедура UMAT; SIMULIA Abaqus; ползучесть; сплошность; поврежденность; асимптотика.

Цитирование. Чаплий Д.В., Степанова Л.В., Белова О.Н. Воздействие аккумуляции повреждений на асимптотическое поведение напряжений в окрестности вершины трещины // Вестник Самарского

Исследование поддержано Российским научном фондом, проект № 21-11-00346.

университета. Естественнонаучная серия. 2023. Т. 29, № 1. С. 47-63. Б01: http://doi.org/10.18287/ 2541-7525-2023-29-1-47-63.

Информация о конфликте интересов: авторы и рецензенты заявляют об отсутствии конфликта интересов.

© Чаплий Д.В., 2023

Чаплий Дмитрий Викторович — аспирант кафедры математического моделирования в механике, Самарский национальный исследовательский университет имени академика С.П. Королева, 443086, Российская Федерация, г. Самара, Московское шоссе, 34.

© Степанова Л.В., 2023

Степанова Лариса Валентиновна — заведующий кафедрой математического моделирования в механике, Самарский национальный исследовательский университет имени академика С.П. Королева, 443086, Российская Федерация, г. Самара, Московское шоссе, 34.

© Белова О.Н., 2023

Белова Оксана Николаевна — ассистент кафедры математического моделирования в механике, Самарский национальный исследовательский университет имени академика С.П. Королева, 443086, Российская Федерация, г. Самара, Московское шоссе, 34.

1. Предварительные сведения

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

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

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

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

Например, в [2] для оценки структурной целостности конструкции при повышенной температуре в реакторе на расплаве соли тория выполнен расчет повреждений при ползучести. Авторы отмечают, что модель, базирующаяся на теории повреждений при ползучести, и метод численного моделирования до сих пор не были предложены для ключевых материалов данного класса реакторов. В [2] исследуются

характеристики повреждений при ползучести сплава UNS N10003 с использованием закона ползучести Нортона и модели повреждений Качанова — Работнова в условиях ползучести. Во-первых, были приняты экспериментальные данные о ползучести сплава UNS N10003 при 650 °C, чтобы соответствовать материальным константам двух моделей. Во-вторых, было проанализировано и обсуждено поведение сплава UNS N10003 при повреждении при ползучести в одно- и многоосном напряженных состояниях. Результаты показали, что модель повреждения в условиях ползучести Качанова — Работнова больше подходит для сплава UNS N10003, чем модель, основанная на классических уравнениях Бейли — Нортона. Наконец, был разработан метод численного моделирования с помощью пользовательской подпрограммы UMAT и впоследствии проверен с помощью анализа МКЭ. Авторы показывают, что результаты конечно-элементного моделирования соответствовали имеющимся теоретическим оценкам. Авторы наглядно демонстрируют, что развитый ими подход представляет собой эффективный метод анализа повреждений в конструкциях в режиме ползучести при повышенной температуре в данном типе реакторов.

Принципы разработки пользовательских моделей материалов UMAT обсуждаются в [3]. Авторы подчеркивают важность основополагающих принципов, поскольку не все известные конституциональные модели материалов были смоделированы и подтверждены в качестве встроенных моделей материалов в рамках обычных МКЭ-решателей. Там, где обнаруживаются новые материалы или даже требуются модификации и усовершенствования существующих встроенных моделей материалов, пользователь МКЭ-пакетов обычно прибегает к вычислительному описанию своих версий конституциональных выражений, чтобы отразить такие изменения.

В [4] авторы представили прогнозирующее моделирование развития повреждений при ползучести материалов путем реализации модифицированных определяющих соотношений Робинсона — Русселье и расширенного метода конечных элементов (XFEM) для решения проблемы разрушения при ползучести в задаче о росте трещины. Разработана новая модель, называемая модифицированной моделью Робинсона — Русселье, для прогнозирования поведения повреждений при ползучести с точки зрения микромеханических повреждений вследствие роста пор в условиях ползучести. Интерфейс модели выполняет неявную схему интеграции в подпрограмме UMAT модуля Abaqus/Standard. Метод радиального возврата используется для интегрирования определяющего уравнения вязкопластичности в конечно-элементной формулировке. Численные модели в 2D и 3D элементах реализуются для определения корректности разработанных подпрограмм, а результаты сравниваются с точным решением для верификации. Кроме того, испытания на ползучесть при растяжении на образце из гладких прутков моделируются и испытываются при постоянной температуре 625 °C с различными уровнями напряжений. Результаты показывают, что максимальные значения напряжений, деформации ползучести и повреждения обнаруживаются вблизи центра образца при растяжении, где наблюдается шейкообразование. Кроме того, результаты сравниваются с известными литературными данными, чтобы проверить и оценить разработанную модель и показать разумное соответствие между обоими результатами. Затем проведенный анализ расширяется путем введения данных о развитии трещин в образце на основе метода XFEM. В результате предложена новая модель, получившая название модифицированной модели Робинсона — Русселье, и результаты сравниваются с результатами, найденными в литературе, которые показали эволюцию роста пустот на пути распространения трещины. Таким образом, доказано [4], что модельное решение обладает потенциалом для прогнозирования поведения повреждений при ползучести с точки зрения роста трещин путем образования, роста и слияния пор на линии продолжения дефекта в конструкциях из пластичных материалов.

В [5] выполнены классический анализ ползучести и моделирование повреждений в непрерывном режиме ползучести для сварного соединения паропровода с исходным материалом 0,5Cr0,5Mo0,25B, металлом сварного шва 2,25CrMo и рассмотрена зона термического воздействия. Оценка срока службы при ползучести с использованием классического анализа деформации ползучести, включающего закон линейного накопления повреждений, часто в значительной степени отличается от практики. Поэтому в исследовании [5] выполнены три вида различного анализа, предполагающих классический анализ ползучести (без деградации материала из-за повреждений), непрерывное повреждение при ползучести при однозонном сварном пересечении и непрерывное повреждение при ползучести при многозонном сварном пересечении. Необходимые формулы для анализа континуальных повреждений при ползучести реализованы в пользовательской подпрограмме материалов (UMAT) и связаны с программным обеспечением ANSYS для конечных элементов. Представлено и обсуждено прогнозируемое распределение эквивалентных напряжений, эволюция повреждений и срок службы зарождения трещин на основе упомянутого выше моделирования. Показано, что прогнозируемый срок службы на основе классического анализа ползучести и с учетом многозонного моделирования значительно превышает срок службы, полученный на основе моделирования повреждений при ползучести. Однако прогнозируемые результаты анализа повреждений при ползучести по одно- и многозонным моделям сопоставимы с разницей примерно в 7,5 %.

Таким образом, к настоящему времени сформировалось устойчивое суждение о необходимости инкорпорирования конституциональных моделей материалов, включающих разнообразные меры поврежденно-сти [1—5]. Тем не менее, несмотря на довольно сильную и достаточную изученность проблемы, которой занимается множество научных учреждений, многие вопросы остаются открытыми и требующими несомненного рассмотрения.

Задача моделирования напряженно-деформированного состояния (НДС) в условиях ползучести являлась фундаментальной проблемой нелинейной механики разрушения и остается актуальной в настоящее время [1-27]. Сегодня особый интерес вызывают исследование разрушения материалов и моделирование поврежденности тел с трещинами в условиях ползучести [6-17].

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

До настоящего времени предложено множество определяющих соотношений и кинетических уравнений для описания эволюции повреждений как совершенно новых [14-16], так и базирующихся на модели поврежденности Качанова — Работнова [23].

Нелинейная механика разрушения была развита с применением 7-интеграла. Впервые введенный Райсом в 1968 году 7-интеграл определяется следующим образом [24]:

7=Г {Шп1 - ац пк £0 (1)

В данном соотношении п — это нормаль к контуру интегрирования Г, Ш — плотность энергии деформации, и* — перемещение.

Известно классическое решение Хатчинсона — Райса — Розенгрена [24; 25], связывающее 7-интеграл и поля НДС в вершине трещины при упругопластическом деформировании. Полученное решение описывается следующими соотношениями:

1

(г, в) = (аЬ) "+1 ^(в, и),

¿j (т,в) = a J n+1 Uij (2)

Ui(r,e) = ai arnf) n+1 r 1+n Uij (9,n).

Эти асимптотические формулы часто называют ХРР-асимптотиками. Пользуясь аналогией Хоффа, распределения напряжений, деформаций и перемещений для случая установившейся ползучести представляются аналогично, но в терминах С * и вычисляются согласно следующим соотношениям:

(г,в) = (ав) "+1 ч (в, и), , ) +

^(г, в) = а(ОЬ) Ь+1 ^(в, и), (3)

, ) + 1 щ(г,в) = аЬ) и*!(в, и)-

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

С* = !с ^*(1'Х2 - ^и* ' (4)

где

W * = оц ¿¿ц. (5)

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

В данном исследовании показано, что на определенном расстоянии независимо от формы вершины трещины выполняется ХРР-решение. Так актуальной и центральной становится задача извлечения

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

Целью настоящего исследования является оценка влияния процесса накопления повреждений на асимптотическое поведение механических характеристик у вершины трещины средствами метода конечных элементов. Для анализа асимптотики полей и е^ в вычислительную процедуру метода конечных элементов были с помощью иМАТ [28; 29] внедрены конституциональные модели материалов, учитывающие процесс накопления повреждений.

2. Моделирование полей, ассоциированных с вершиной трещины в режиме ползучести

В пакете моделирования 81МИЫА АЪадив, осуществляющем расчеты методом конечных элементов, проведено моделирование плоской пластины с центральной трещиной в условиях ползучести. Пластина имеет размеры 100 х 100 мм, длина трещины равна 10 мм. В первой серии расчетов трещина моделируется острой. Пластина находится в условиях одноосного растяжения.

Для моделирования ползучести избран степенной закон Бейли — Нортона, имеющий следующий вид:

3 ! /З

éij = 2 Bae- Sij ' = y 2 Sij Sj ' (6)

где èij — компоненты скорости деформации ползучести, а^ — компоненты напряжений Коши, ае — интенсивность напряжения, n,B — константы материала; Sij = aij — аииôij/3 — компоненты девиатора напряжений, где ôij — символ Кронекера. Это определяющее соотношение было внедрено в расчет с помощью пользовательской подпрограммы UMAT.

Материальные константы n,B в соотношении (6) — это параметры материала, в общем случае зависящие от температуры. В качестве материала пластины выбрана сталь, имеющая следующие упругие модули: модуль Юнга E = 210000 Н/мм2, коэффициент Пуассона v равен 0.3. Проведен ряд расчетов с различными материальными константами B и n. Принятые значения констант и условия экспериментальных исследований приведены в табл. 1.

Таблица 1

Условия проведения численных расчетов

Table 1

Conditions for making settlement calculations

№ Модуль Коэффициент B, n Нагрузка

модели Юнга Пуассона (Н/мм2)-п(ч)-1 P, Н/мм2

1.1 210000 0.3 1 • 10-10 2 20

1.2 210000 0.3 1 • 10-12 3 20

1.3 210000 0.3 1 • 10-15 4 20

1.4 210000 0.3 1 • 10-16 5 20

Обычно показатель ползучести п принимает значения от 3 до 8 для металлов и сплавов. Для отдельных современных высокопрочных сплавов достигает 18 [15; 22; 23]. Изображение модели пластины со схемой приложенной нагрузки представлено на рис. 1.

Сетка построена с использованием сингулярных конечных элементов в окрестности вершины трещины. Концентрические окружности поделены на 72 сектора, таким образом, раствор секторального элемента равен 5°. Тип элементов сетки был избран СРБ4. Полученная конечно-элементная сетка представлена на рис. 2.

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

Распределение компоненты тензора напряжений приведено вдоль радиальной траектории, построенной под углом в = 0, показанной на рис. 3, и построено в двойных логарифмических координатах.

Рис. 1. Изображение модели со схемой приложенной нагрузки Fig. 1. Scheme of the model geometry and applied loads

Рис. 2. Конечно-элементная сетка в модели Fig. 2. Typical finite element mesh of the model

Рис. 3. Путь, вдоль которого строятся радиальные распределения напряжений Fig. 3. The path along which stress radial distributions are constructed

На рис. 4 изображены итоги вычислений в системе 81МиЫЛ ЛЪадив для четырех различных моделей с острой центральной трещиной со свойствами в соответствии с табл. 1.

Поля компонент тензора напряжений <22 приведены в различные моменты времени. Общее время моделирования равно 35 000 часов. Различными оттенками отмечены радиальные распределения поля напряжений в разные моменты времени.

d

Рис. 4. Результаты конечно-элементного анализа: компоненты тензора напряжений а22 в условиях ползучести c течением времени вдоль радиальной траектории (в = 0). Результаты для: a — модели 1.1,

b — модели 1.2, c — модели 1.3, d — модели 1.4 Fig. 4. Results of FEM: the stress tensor components a22 under creep conditions along the radial trajectory (в = 0). Results for: a — model 1.1, b — model 1.2, c — model 1.3, d — model 1.4

a

b

с

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

Для каждого графика на рис. 4 построена прямая черная линия, аппроксимирующая значения компоненты тензора напряжений <22 в последний момент времени в зоне ползучести. Полученные значения коэффициента наклона этих прямых соответствуют значениям, полученным из соотношения к = .

Линия, аппроксимирующая значения компоненты тензора напряжений <22 в зоне упругости, имеет при этом коэффициент наклона равный к = 2, что также соответствует аналитическому решению. Это показывает, что асимптотическое решение ХРР присутствует в расчетах для острой трещины.

Многочисленные исследования показали, что расчеты напряженно-деформированного состояния в условиях ползучести необходимо проводить в образцах с трещинами, имеющими закругленные вершины [12; 13; 27]. Ввиду этого следующая серия вычислений проведена для трещин с закругленными концами. Трещина моделируется в виде разреза с закругленными вершинами, радиус закругления составляет 0.001 мм. В настоящем исследовании проведен ряд расчетов с различными материальными константами В и п. Принятые значения констант и условия экспериментальных исследований приведены в табл. 2.

Таблица 2

Условия проведения численных расчетов

Table 2

Conditions for making settlement calculations

№ модели Модуль Юнга Коэффициент Пуассона B, (Н/мм2)-п(ч)-1 n Нагрузка P, Н/мм2

2.1 210000 0.3 1 • 10-10 2 20

2.2 210000 0.3 1 • 10-12 3 20

2.3 210000 0.3 1 • 10-15 4 20

Сетка в окрестности вершины трещины построена с использованием сингулярных конечных элементов. Концентрические окружности поделены на три области: перед вершиной дуга окружности величиной 180° разбита на 36 секторов, таким образом, раствор секторального элемента равен 5°, две другие части поделены на 17 секторов. Тип элементов сетки был избран CPS4. Полученная конечно-элементная сетка представлена на рис. 5.

Рис. 5. Конечно-элементная сетка в модели с закругленной трещиной Fig. 5. Typical finite element mesh of the model

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

На рис. 7 изображены поля компонент тензора напряжений а22 в логарифмических координатах с течением времени вдоль траектории, показанной на рис. 6. Различными оттенками отмечены радиальные распределения поля напряжений в разные моменты времени.

Для каждого графика построена черная линия, аппроксимирующая значения компоненты тензора напряжений а22 в последний момент времени в зоне ползучести. Полученные значения коэффициента наклона этих прямых находятся в согласии со значениями, соответствующим формуле к = n+j. Можно отметить, что в области, близкой к вершине закругленной трещины, асимптотика меняется, в отличие от предыдущих расчетов с острой трещиной. Это связано с тем, что такая конфигурация трещины приводит к изменению сингулярности вблизи трещины.

Рис. 6. Путь, вдоль которого строятся радиальные распределения напряжений Fig. 6. The path along which stress radial distributions are constructed

Линия, аппроксимирующая значения компоненты тензора напряжений <22 в зоне упругости, имеет при этом коэффициент наклона равный к = 1, что также соответствует аналитическому решению.

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

3. Анализ влияния прогрессивного накопления

микроповреждений у кончика трещины на асимптотику

Задачей последующего исследования выступает анализ влияния развитой эволюции поврежденности вблизи вершины закругленной трещины в предположении реализации плоского напряженного состояния для степенного закона установившейся ползучести и кинетического уравнения Качанова — Работнова на асимптотику напряжений.

Степенной закон ползучести, следуя классической процедуре [20], был дополнен кинетическим уравнением накопления поврежденности Качанова — Работнова и математическая модель представляется в следующем виде: )

-=3 Ч1Г • *=—^^Т- <7>

где ф — параметр сплошности Качанова [18]; а^к — гидростатическое напряжение, — максимальное главное напряжение, аедЬ = аае + в<1 + (1 — а — в)акк — эквивалентное напряжение, а константы а и в находятся экспериментально. Неповрежденному, целостному материалу соответствует значение параметра ф = 1, а ф = 0 означает, что материал полностью исчерпал несущую способность. Вместо параметра сплошности часто используется параметр поврежденности Работнова ш [20], связанный с параметром сплошности выражением ш = 1 — ф.

Для конечно-элементного расчета необходимо задать постоянные материала А и т, фигурирующие в эволюционном уравнении 7. Обычно т выбирают таким образом, чтобы выполнялось соотношение т « 0.7п, полученное эмпирическим путем [22].

В ходе исследования проведен ряд расчетов с различными материальными константами А,В,п и т. Принятые значения констант и условия экспериментальных исследований приведены в табл. 3.

Постоянные материала в эволюционном уравнении и определяющем соотношении (7) выбирались на основании экспериментальных данных, описанных в работах [15; 23], где приведены значения материальных констант для большого ряда металлов и сплавов. С использованием подпрограммы, являющейся процедурой ИМАТ комплекса 81МИЫА АЪадив, конституциональные соотношения степенного закона ползучести и эволюционное уравнение (7) были внесены в вычислительную схему метода конечных элементов пакета Abaqus/Standart. Это дало возможность найти распределения параметра сплошности и механических полей.

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

b

c

Рис. 7. Результаты конечно-элементного анализа: компоненты тензора напряжений а22 в условиях ползучести c течением времени. Результаты приведены для: a — модели 2.1, b — модели 2.2, c — модели 2.3 Fig. 7. Results of FEM: the stress tensor components a22 under creep conditions. Results for: a — model 2.1,

b — model 2.2, c — model 2.3

Таблица 3

Условия проведения численных расчетов

Table 3

Conditions for making settlement calculations

№ модели Модуль Юнга Коэффициент Пуассона B, (Н/мм2)-п(ч)-1 n A, (Н/мм2 )-п(ч)-1 m Нагрузка P, Н/мм2

3.1 210000 0.3 1 • 10-13 3 1•10-1 2.5 10

3.2 210000 0.3 1 • 10-16 4 1•10-1 2.8 10

3.3 210000 0.3 1 • 10-19 5 1•10-1 3.5 10

На рис. 8 изображены поля компоненты тензора напряжений <22 в двойных логарифмических координатах с течением времени вдоль траектории, показанной на рис. 6 и построенной по прямой под

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

Рис. 8. Результаты конечно-элементного анализа: компоненты тензора напряжений а22 в условиях ползучести при действии процесса накопления поврежденности c течением времени. Результаты для:

a — модели 3.1, b — модели 3.2, c — модели 3.3 Fig. 8. Results of FEM: the stress tensor components a22 under creep condition under the action of damage accumulation over time. Results for: a — model 3.1, b — model 3.2, c — model 3.3

с

Для каждого графика построена регрессионная прямая, отмеченная на графиках черной линией, аппроксимирующая значения компоненты тензора напряжений <22 в последний рассматриваемый момент времени в зоне ползучести. Полученные значения коэффициента наклона этой прямой существенно отличаются от значений, соответствующих формуле k = Пру, а значит и асимптотике решения ХРР. Здесь также можно отметить, что в области, близкой к вершине закругленной трещины, асимптотика не наблюдается, в отличие от расчетов с острой трещиной. Результаты на рис. 8 отличаются от распределений, полученных на расчетах при действии ползучести без учета накопления поврежденности, как на рис. 7. По графикам на рис. 8 можно заметить, каким образом процесс накопления поврежденности изменяет асимптотику напряжений в окрестности вершины трещины.

На рис. 9 показано распределения сплошности с течением времени: a-c — состояние в 922.9, 3 923, 7 923 часа, d-f — состояние в 222.9, 9 923, 3 5000 часов.

На рис. 10 представлены картины распределения сплошности вдоль траектории, показанной на рис. 6, с течением времени. Различными оттенками отмечены радиальные распределения сплошности в разные моменты времени.

d

Рис. 9. Распределение сплошности c течением времени в a-c модели 3.1, d-f в модели 3.3 Fig. 9. Distribution of continuity over time a-c for model 3.1, d-f for model 3.3

a

Рис. 10. Результаты конечно-элементного анализа: распределение сплошности в условиях ползучести с учетом накопления поврежденности c течением времени. Результаты для: a — модели 3.1, b — модели 3.2,

c — модели 3.3

Fig. 10. Results of FEM: distribution of continuity under creep conditions, taking into account the accumulation of damage over time. Results for: a — model 3.1, b — model 3.2, c — model 3.3

Из рис. 9, 10 можно заметить, каким образом для разных значений материальных констант меняется поле сплошности.

Заключение

В ходе исследования выполнено конечно-элементное моделирование одноосного растяжения двумерной пластины с центральной трещиной, находящейся в условиях ползучести с учетом поврежденности в связанной постановке. Моделирование выполнено в комплексе 81МиЫЛ ЛЪадив с применением пользовательской процедуры иМЛТ. Степенной закон ползучести с помощью пользовательской процедуры был дополнен кинетическим уравнением накопления поврежденности Качанова — Работнова в связанной постановке. В результате моделирования получены распределения напряжений, деформаций и сплошности в условиях ползучести с учетом накопления поврежденности с течением времени. Приведены асимптотики распределения компоненты тензора напряжений < 22 для моделей как с острой трещиной, так и с вырезом с закругленными вершинами. Моделирование выполнено для различных значений материальных констант. Приведено радиальное распределение сплошности вдоль траектории, построенной под углом в = 0 к горизонтальной оси х. Показано, что асимптотика в упругой области и в области развитой ползучести в первых двух сериях расчетов соответствует аналитическим решениям, однако в случае учета накопления поврежденности эта асимптотика не наблюдается.

Литература

[1] Hayhurst D.R. Materials Data Bases and Mechanisms-Based Constitutive Equations for Use in Design // In: Altenbach, H., Skrzypek, J.J. (Eds.) Creep and Damage in Materials and Structures. International Centre for Mechanical Sciences, vol 399. Vienna: Springer-Verlag, 1999, pp. 167-208. DOI: http://doi.org/10.1007/978-3-7091-2506-9_4.

[2] Wang X.-Y., Zhang X.-C., Zhu S.-F. Creep damage characterization of UNS N10003 alloys based on a numerical simulation using the Norton creep law and Kachanov-Rabotnov creep damage model // Nuclear Science and Techniques. 2019. Vol. 30. Article number 65. DOI: http://dx.doi.org/10.1007/s41365-019-0586-2.

[3] Okereke M., Keates S. Material Response: Constitutive Models and Their Implementation. // In: Finite Element Applications. Springer Tracts in Mechanical Engineering. Cham: Springer, 2018, pp. 363-436. DOI: http://doi.org/10.1007/978-3-319-67125-3_10.

[4] Ahmad M.I.M, Akbar M., Abdullah N.A. Development of constitutive creep damage-based modified Robinson-Rousselier (MRR) model with XFEM for void-crack relation in ductile materials // Mechanics of Time-Dependent Materials. 2022. DOI: http://doi.org/10.1007/s11043-022-09540-5.

[5] Hosseini-Toudeshky H., Jannnesari M. An investigation on creep life assessment of welded steam pipeline intersection using classical and progressive damage analyses // Welding in the World, 2022, Vol. 66, pp. 1653-1664. DOI: http://doi.org/10.1007/s40194-022-01324-2.

[6] Степанова Л.В. Компьютерное моделирование процессов накопления повреждений в твердых телах с трещинами с помощью пользовательской процедуры UMAT вычислительного комплекса Simulia Abaqus // Вестник Пермского национального исследовательского политехнического университета. Сер.: Механика. 2018. №. 3. С. 71-86. DOI: http://doi.org/10.15593/perm.mech/2018.3.08. EDN: https://www.elibrary.ru/yljexz.

[7] Turkova V.A., Stepanova L.V. Evaluation of damage accumulation zone in the vicinity of the crack tip: FEM analysis via UMAT procedure // Journal of Physics: Conference Series. 2018. Vol. 1096. P. 012157. DOI: http://doi.org/10.1088/1742-6596/1096/1/012157.

[8] Шлянников В.Н., Туманов А.В. Силовая и деформационная модели поврежденности и разрушения при ползучести // Физическая мезомеханика. 2018. Т. 21, № 3. C. 70-85. DOI: http://doi.org/10.24411/1683-805X-2018-13008. EDN: https://www.elibrary.ru/xrgsgd.

[9] Shlyannikov V.N. Crack tip fields and fracture resistance parameters based on strain gradient plasticity // International Journal of Solids and Structures. 2021. Vols. 208-209, pp. 63-82. DOI: https://doi.org/10.1016/j.ijsolstr.2020.10.015.

[10] Шлянников В.Н. Решение задач нелинейного деформирования и разрушения материалов при сложном напряженном состоянии // Физическая мезомеханика. 2012. Т. 15, № 1. C. 57-67. URL: https://www.elibrary.ru/item.asp?id=17428838. EDN: https://www.elibrary.ru/orkhbf.

Shlyannikov V. N., Tumanov A. V. Creep damage and stress intensity factor assessment for plane multi-axial and three-dimensional problems // International Journal of Solids and Structures. 2018. Vol. 150. Pp. 166-183. DOI: https://doi.org/10.1016/j.ijsolstr.2018.06.009. EDN: https://www.elibrary.ru/xwikmx.

McMeeking R.M. Finite deformation analysis of crack-tip opening in elastic-plastic materials and implications for fracture // Journal of the Mechancis and Physics of Solids. 1977. Vol. 25, № 5. Pp. 357-381. DOI: http://doi.org/10.1016/0022-5096(77)90003-5.

Shlyannikov V.N., Kislova S.Y. Parameters of deformation mixed modes with account of crack tip curvature // Strength of Materials. 2010. Vol. 42, № 6. Pp. 660-674. DOI: http://doi.org/10.1007/s11223-010-9254-9.

Meng Li., Chen W., Yan Y., Kitamura T., Feng. M. Modelling of creep and plasticity deformation considering creep damage and kinematic hardening // Engineering Fracture Mechanics. 2019. Vol. 218. P. 106582. DOI: http://dx.doi.org/10.1016/j.engfracmech.2019.106582.

Meng Q., Zhenqing H. Creep damage models and their applications for crack growth analysis in pipes: A review // Engineering Fracture Mechanics. 2019. Vol. 205. Pp. 547-576. DOI: http://dx.doi.org/10.1016/j.engfracmech.2015.09.055.

Murakami S. Continuum Damage Mechanics. A Continuum Mechanics Approach to the Analysis of Damage and Fracture. Springer Dordrecht, 2012. 423 p. Available at: https://doi.org/10.1007/978-94-007-2666-6.

Wen J.F., Tu S.T., Gao X.L., Reddy J.N. Simulations of creep crack growth in 316 stainless steel using a novel creep-damage model // Engineering Fracture Mechanics. 2012. Vol. 98. Pp. 169-184. DOI: http://doi.org/10.1016/j.engfracmech.2012.12.014.

Kachanov L.M. Force and deformation models of damage and fracture during creep // Physical Mesomechanics. 2018. Vol. 21, № 3. Pp. 70-85.

Качанов Л.М. Теория ползучести. Mосква: Гос. изд-во физико-математической литературы, 1960. 455 с. URL: https://lib-bkm.ru/13835.

Работнов Ю.Н. Ползучесть элементов конструкций. Москва: Наука, 2014. 752 c. URL: https://lib-bkm.ru/13795.

Локощенко А.М., Фомин Л.В., Теруад В.В., Басалов Ю.Г., Агабабян В.С. Ползучесть и длительная прочность металлов при нестационарных сложных напряженных состояниях (обзор) // Вестник Самарского государственного технического университета. Сер.: Физико-математические науки. 2020. Т. 24, № 2. C. 275-318. DOI: http://doi.org/10.14498/vsgtu1765. EDN: https://www.elibrary.ru/oqccvc.

Бойл Дж., Спенс Дж. Анализ напряжений в конструкциях при ползучести. Москва: Мир. 1986. 360 с.

Riedel H. Fracture at High Temperatures. Berlin, Heidelberg: Springer-Verlag, 1987. 418 p. DOI: http://doi.org/10.1007/978-3-642-82961-1.

Rice J.R., Rosengren G.F. Plane strain deformation near a crack tip in a power-law hardening material // Journal of The Mechanics and Physics of Solids. 1968. Vol. 16, № 1. Pp. 1-12. URL: http://esag.harvard.edu/rice/016_RiceRosengren_CrackSing_JMPS68.pdf.

Hutchinson J.W. Singular behaviour at the end of a tensile crack in a hardening material // Journal of The Mechanics and Physics of Solids. 1968. Vol. 16, № 1. Pp. 13-31. URL: https://groups.seas.harvard.edu/hutchinson/papers/312.pdf.

Hutchinson J.W. Plastic stress and strain fields at a crack tip // Journal of The Mechanics and Physics of Solids. 1968. Vol. 16, № 5. Pp. 337-342. URL: https://groups.seas.harvard.edu/hutchinson/papers/313.pdf.

Naumenko K., Altenbach H. Modelling of Creep for Structural Analysis. Berlin, Heidelberg: Springer. 2007. 220 p. DOI: http://dx.doi.org/10.1007/978-3-540-70839-1_1.

USER MATERIAL IN ABAQUS. URL: https://abaqus-docs.mit.edu/2017/English/SIMACAESUBRefMap/simasub-c-umat.htm.

Lecture 6. Writing a UMAT or VUMAT. URL: https://imechanica.org/files/Writing%20a%20UMAT.pdf.

DOI: 10.18287/2541-7525-2023-29-1-47-63 Submited: 03.03.2023

Revised: 14.04.2023 Accepted: 30.05.2023

D.V. Chapliy

Samara National Research University, Samara, Russian Federation E-mail: [email protected]. ORCID: https://orcid.org/0000-0001-9510-3659

L.V. Stepanova

Samara National Research University, Samara, Russian Federation E-mail: [email protected]. ORCID: https://orcid.org/0000-0002-6693-3132

O.N. Belova

Samara National Research University, Samara, Russian Federation E-mail: [email protected]. ORCID: https://orcid.org/0000-0002-4492-223X

EFFECT OF DAMAGE ACCUMULATION ON THE ASYMPTOTIC BEHAVIOR OF STRESSES AHEAD THE CRACK TIP2

ABSTRACT

The subject of this study is the analysis of mechanical fields associated with a crack tip under creep conditions, taking into account the phenomenon of damage accumulation. The objective of the study is to perform finite element modeling, using the SIMULIA Abaqus software package, of uniaxial tension of a plate with a central horizontal crack under creep conditions, taking into account damage accumulation. For numerical simulation of creep, the Bailey-Norton power law is used. The power law of creep with the help of the user procedure UMAT (User Material) of the SIMULIA Abaqus package was supplemented with the Kachanov-Rabotnov kinetic equation of damage accumulation in a related formulation. In the calculation scheme of finite elements, the crack tip was modeled as a mathematical notch and as a notch with a finite radius of curvature. As a result of the calculations, the distributions of stresses, strains, and continuity under creep conditions were obtained, taking into account the accumulation of damage over time. Radial distributions of continuity, stresses, and strains are plotted over time at various distances from the crack tip. The subject of the study was the consideration of the asymptotic of the stress distribution. As a result of the study, it is shown that in the elastic region the asymptotic corresponds to the distribution under the elastic regime, and in the creep zone the asymptotics of Hutchinson, Rice and Rosengren (HRR-solution) is satisfied for different exponents n of the power law of creep.

A comparison is made of the radial stress distributions in modeling without taking into account damage and in the case of taking into account damage accumulation. It is shown that the presence of damage significantly changes the asymptotics of the stress field in the vicinity of the crack tip.

Key words: user procedure UMAT; SIMULIA Abaqus; creep; damage; continuity; asymptotic.

Citation. Chapliy D.V., Stepanova L.V., Belova O.N. Effect of damage accumulation on the asymptotic behavior of stresses ahead the crack tip. Vestnik Samarskogo universiteta. Estestvennonauchnaia seriia = Vestnik of Samara University. Natural Science Series, 2023, vol. 29, no. 1, pp. 47-63. DOI: http://doi.org/10.18287/2541-7525-2023-29-1-47-63. (In Russ.)

Information about the conflict of interests: authors and reviewers declare no conflict of interests.

© Chapliy D.V., 2023

Dmitriy V. Chapliy — postgraduate student of the Department of Mathematical Modeling in Mechanics, Samara National Research University, 34, Moskovskoye shosse, Samara, 443086, Russian Federation.

©c Stepanova L.V., 2023

Larisa V. Stepanova — Doctor of Physical and Mathematical Sciences, associate professor, head of the Department of Mathematical Modeling in Mechanics, Samara National Research University, 34, Moskovskoye shosse, 443086, Russian Federation.

©c Belova O.N., 2023

Oksana N. Belova — assistant lecturer of the Department of Mathematical Modeling in Mechanics, Samara National Research University, 34, Moskovskoye shosse, Samara, 443086, Russian Federation.

2The research is supported by the Russian Scientific Fund, project № 21-11-00346.

References

[1] Hayhurst D.R. Materials Data Bases and Mechanisms-Based Constitutive Equations for Use in Design. In: Altenbach, H., Skrzypek, J.J. (Eds.) Creep and Damage in Materials and Structures. International Centre for Mechanical Sciences, vol 399. Vienna: Springer-Verlag, 1999, pp. 167-208. DOI: http://doi.org/10.1007/978-3-7091-2506-9_4.

[2] Wang X.-Y., Zhang X.-C., Zhu S.-F. Creep damage characterization of UNS N10003 alloys based on a numerical simulation using the Norton creep law and Kachanov-Rabotnov creep damage model. Nuclear Science and Techniques, 2019, vol. 30. Article number 65. DOI: http://dx.doi.org/10.1007/s41365-019-0586-2.

[3] Okereke M., Keates S. Material Response: Constitutive Models and Their Implementation. In: Finite Element Applications. Springer Tracts in Mechanical Engineering. Cham: Springer, 2018, pp. 363-436. DOI: http://doi.org/10.1007/978-3-319-67125-3_10.

[4] Ahmad M.I.M, Akbar M., Abdullah N.A. Development of constitutive creep damage-based modified Robinson-Rousselier (MRR) model with XFEM for void-crack relation in ductile materials. Mechanics of Time-Dependent Materials, 2022. DOI: http://doi.org/10.1007/s11043-022-09540-5.

[5] Hosseini-Toudeshky H., Jannnesari M. An investigation on creep life assessment of welded steam pipeline intersection using classical and progressive damage analyses. Welding in the World, 2022, vol. 66, pp. 1653-1664. DOI: http://doi.org/10.1007/s40194-022-01324-2.

[6] Stepanova L.V. Computational simulation of the damage accumulation processes in cracked solids by the user procedure UMAT of SIMULIA Abaqus. PNRPU Mechanics Bulletin, 2018, no. 3, pp. 71-86. DOI: http://doi.org/10.15593/perm.mech/2018.3.08. EDN: https://www.elibrary.ru/yljexz. (In Russ.)

[7] Turkova V.A., Stepanova L.V. Evaluation of damage accumulation zone in the vicinity of the crack tip: FEM analysis via UMAT procedure. Journal of Physics: Conference Series, 2018, vol. 1096, p. 012157. DOI: http://doi.org/10.1088/1742-6596/1096/1/012157.

[8] Shlyannikov V.N., Tumanov A.V. Force and deformation models of damage and fracture during creep. Physical Mesomechanics, 2018, vol. 21, № 3, pp. 70-85. DOI: https://doi.org/10.24411/1683-805X-2018-13008. EDN: https://www.elibrary.ru/xrgsgd. (In Russ.)

[9] Shlyannikov V.N. Crack tip fields and fracture resistance parameters based on strain gradient plasticity. International Journal of Solids and Structures, 2021, vols. 208-209, pp. 63-82. DOI: http://doi.org/10.1016/j.ijsolstr.2020.10.015.

[10] Shlyannikov V.N. Solution of nonlinear strain and fracture problems of materials in complex stress states. Physical Mesomechanics, 2012, vol. 15, no. 1, pp. 57—67. Available at: https://www.elibrary.ru/item.asp?id=17428838. EDN: https://www.elibrary.ru/orkhbf. (In Russ.)

[11] Shlyannikov V.N., Tumanov A.V. Creep damage and stress intensity factor assessment for plane multi-axial and three-dimensional problems. International Journal of Solids and Structures, 2018, vol. 150, pp. 166-183. DOI: http://doi.org/10.1016/j.ijsolstr.2018.06.009. EDN: https://www.elibrary.ru/xwikmx.

[12] McMeeking R.M. Finite deformation analysis of crack-tip opening in elastic-plastic materials and implications for fracture. Journal of the Mechancis and Physics of Solids, 1977, vol. 25, № 5, pp. 357-381. DOI: http://doi.org/10.1016/0022-5096(77)90003-5.

[13] Shlyannikov V.N., Kislova S.Y. Parameters of deformation mixed modes with account of crack tip curvature. Strength of Materials, 2010, vol. 42, № 6, pp. 660-674. DOI: http://doi.org/10.1007/s11223-010-9254-9.

[14] Meng Li., Chen W., Yan Y., Kitamura T., Feng. M. Modelling of creep and plasticity deformation considering creep damage and kinematic hardening. Engineering Fracture Mechanics, 2019, vol. 218, p. 106582. DOI: http://dx.doi.org/10.1016/j.engfracmech.2019.106582.

[15] Meng Q., Zhenqing H. Creep damage models and their applications for crack growth analysis in pipes: A review. Engineering Fracture Mechanics, 2019, vol. 205, pp. 547-576. DOI: http://dx.doi.org/10.1016/j.engfracmech.2015.09.055.

[16] Murakami S. Continuum Damage Mechanics. A Continuum Mechanics Approach to the Analysis of Damage and Fracture. Springer Dordrecht, 2012. 423 p. DOI: https://doi.org/10.1007/978-94-007-2666-6.

[17] Wen J.F., Tu S.T., Gao X.L., Reddy J.N. Simulations of creep crack growth in 316 stainless steel using a novel creep-damage model. Engineering Fracture Mechanics, 2012, vol. 98, pp. 169-184. DOI: http://doi.org/10.1016/j.engfracmech.2012.12.014.

[18] Kachanov L.M. Force and deformation models of damage and fracture during creep. Physical Mesomechanics, 2018, vol. 21, № 3, pp. 70-85.

[19] Kachanov L.M. Creep theory. Moscow: Gosudarstvennoe izdatel'stvo fiziko-matematicheskoi literatury, 1960, 455 p. Available at: https://lib-bkm.ru/13835. (In Russ.)

[20] Rabotnov Yu.N. Creep problems in structural members. Moscow: Nauka, 2014, 752 p. Available at: https://lib-bkm.ru/13795. (in Russ.)

[21] Lokoshchenko A.V., Fomin L.V., Teraud W.V., Basalov Y.G., Agababyan V.S. Creep and long-term strength of metals under unsteady complex stress states (review). Journal of Samara State Technical University. Ser. Physical and Mathematical Sciences, 2020, vol. 24, №2, pp. 275-318. DOI: http://doi.org/10.14498/vsgtu1765. EDN: https://www.elibrary.ru/oqccvc. (In Russ.)

[22] Boyle J.T, Spence J. Stress analysis for creep. Moscow: Mir, 1986, 360 p. (In Russ.)

[23] Riedel H. Fracture at High Temperatures. Berlin, Heidelberg: Springer-Verlag, 1987. 418 p. DOI: http://doi.org/10.1007/978-3-642-82961-1.

[24] Rice J.R., Rosengren G.F. Plane strain deformation near a crack tip in a power-law hardening material. Journal of The Mechanics and Physics of Solids, 1968, vol. 16, № 1, pp. 1-12. Available at: http://esag.harvard.edu/rice/016_RiceRosengren_CrackSing_JMPS68.pdf.

[25] Hutchinson J.W. Singular behaviour at the end of a tensile crack in a hardening material. Journal of The Mechanics and Physics of Solids, 1968, vol. 16, № 1, pp. 13-31. Available at: https://groups.seas.harvard.edu/hutchinson/papers/312.pdf.

[26] Hutchinson J.W. Plastic stress and strain fields at a crack tip. Journal of The Mechanics and Physics of Solids, 1968, vol. 16, № 5, pp. 337-342. Available at: https://groups.seas.harvard.edu/hutchinson/papers/313.pdf.

[27] Naumenko K., Altenbach H. Modelling of Creep for Structural Analysis. Berlin, Heidelberg: Springer, 2007. 220 p. DOI: http://dx.doi.org/10.1007/978-3-540-70839-1_1.

[28] USER MATERIAL IN ABAQUS. Available at: https://abaqus-docs.mit.edu/2017/English/SIMACAESUBRefMap/ simasub-c-umat.htm.

[29] Lecture 6. Writing a UMAT or VUMAT. Available at: https://imechanica.org/files/Writing%20a%20UMAT.pdf.