Прикладная математика и механика, 2020, T. 84, № 5, стр. 663-674

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

Основные способы их решения для случая многоосного циклического нагружения связаны с проведением испытаний образцов, обобщением закономерностей, установленных для одноосных нагружений и описываемых усталостными S–N кривыми типа Веллера и соотношениями типа Баскина [1].

Большое количество критериев (stress-based criteria) основано на прямом обобщении S–N кривых, построенных по результатам усталостных испытаний, и они делятся на две большие группы. В первую группу входят критерии, которые включают амплитуды инвариантных характеристик напряженного состояния в цикле нагружения – октаэдрических напряжений, главных напряжений и т.п. [2, 3], а во вторую группу входят те, которые учитывают значения амплитуд касательных и/или нормальных напряжений на так называемой критической плоскости [4–9]. Как правило, эта плоскость определяется из условия максимума амплитуд касательных, нормальных напряжений или их определенной комбинации по площадкам всевозможных ориентаций. Обзоры по этой тематике представлены в работах [10–13].

Для того чтобы исследовать процессы развития зон усталостных повреждений также существует два подхода. Первый основан на классических представлениях механики разрушения и связывает условия развития усталостных трещин при увеличении числа циклов с амплитудами коэффициентов интенсивности напряжений в вершине трещины. Основное уравнение было предложено Парисом [14], существует большое количество его модификаций [15, 16].

Второй подход использует представления теории повреждаемости, восходящей к работам [17, 18] и развитый впоследствии в [19, 20]. Он применялся [21–23] в приложении к задачам циклического нагружения и усталостного разрушения.

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

Для выделения различных режимов усталостного разрушения используется схема мультирежимной амплитудной усталостной кривой, представленная на рис. 1. Вплоть до значения N ~ 10³ реализуется режим повторно-статического нагружения с амплитудой, мало отличающейся от статического предела прочности ${{\sigma }_{B}}$. Далее левая часть бимодальной усталостной кривой (кривая Веллера) описывает режимы МЦУ–МНЦУ вплоть до N ~ 10⁷ и значений амплитуды порядка предела усталости ${{\sigma }_{u}}$. Затем начинается зона смены механизмов разрушения и дальнейшее падение усталостной прочности, начиная с величин N ~ 10⁸, до нового предельного значения ${{\tilde {\sigma }}_{u}}$ в соответствии с правой ветвью бимодальной усталостной S–N кривой. Эта ветвь описывает режим СВМУ [24].

Отметим, что в настоящее время развиваются представления о явном разделении классической ветви Веллера еще на две части – собственно, МЦУ и МНЦУ. Граница этой переходной области определяется не значением величины N, а величиной амплитуды нагружения, равной пределу текучести материала ${{\sigma }_{T}}$ [25], так как при этом также происходит смена физического механизма усталостного разрушения.

Граница повторно-статического диапазона N ~ 10³ довольно условна. Она уточняется в современных исследованиях [25] в зависимости от прочностных и пластических характеристик материала.

Тем не менее, в данной работе мы сохраним изложение предлагаемой модели развития повреждаемости, основанное на вышеописанной схеме.

Для того чтобы согласовать модель с известными критериями многоосного усталостного разрушения, был выбран основанный на напряженном состоянии (stress-based) критерий, который описывает усталостное разрушение, связанное с развитием микротрещин нормального отрыва. Это модифицированный критерий Smith–Watson–Topper (SWT) [26, 27], в котором определяющую роль в развитии усталостных повреждений играют амплитуды максимальных растягивающих напряжений. В исходном изложении этот и подобные ему stress-based критерии уже применялись для описания режимов МЦУ и МНЦУ. Был предложен [24] подход для обобщения многоосных критериев разрушения при описании правых ветвей усталостных кривых в режиме СВМУ, использующий опорные точки каждой из ветвей и обратную степенную зависимость от числа циклов N для выхода на асимптоту предела усталости. В данной работе для обобщения критерия SWT на случай СВМУ был также использован этот аналитический подход.

1. Кинетическое уравнение для повреждаемости в режиме МЦУ–МНЦУ.

1.1. Критерий многоосного усталостного разрушения в режиме МЦУ–МНЦУ с развитием микротрещин нормального отрыва (stress-based SWT), соответствующий левой ветви бимодальной усталостной кривой (рис. 1) имеет вид:

где ${{\sigma }_{1}}$ – наибольшее главное напряжение, $\Delta {{\sigma }_{1}}$ – размах наибольшего главного напряжения за цикл, $\Delta {{\sigma }_{1}}{\text{/}}2$ – его амплитуда. Из условия повторного статического разрушения вплоть до значений N ~ 10³ по методике [13] можно получить, что ${{\sigma }_{L}} = {{10}^{{3{{\beta }_{{LH}}}}}}({{\sigma }_{B}} - {{\sigma }_{u}})$. В соответствии с выбранным критерием к разрушению приводят только растягивающие напряжения, поэтому в него входит величина $\left\langle {{{\sigma }_{{{{1}_{{\max }}}}}}} \right\rangle = {{\sigma }_{{{{1}_{{\max }}}}}}H({{\sigma }_{{{{1}_{{\max }}}}}})$, где ${{\sigma }_{B}}$ – статический предел прочности материала, ${{\sigma }_{u}}$ – классический предел усталости материала при реверсивном цикле (коэффициент асимметрии цикла $R = {{\sigma }_{{\min }}}{\text{/}}{{\sigma }_{{\max }}} = - 1$), ${{\beta }_{{LH}}}$ – степенной показатель левой ветви бимодальной усталостной кривой, $H(x)$ – функция Хевисайда.

Из критерия усталостного разрушения получаем число циклов до разрушения при однородном напряженном состоянии:

1.2. Для описания процесса развития усталостной поврежденности в режиме МЦУ–МНЦУ вводится функция повреждаемости $0 \leqslant \psi (N) \leqslant 1$, которая описывает процесс постепенного циклического (усталостного) разрушения материала. При $\psi = 1$ материальная частица считается полностью разрушенной. Ее модули Ламе становятся равными нулю. Функция повреждаемости $\psi $ в зависимости от числа циклов нагружения $N$ для режима МЦУ-МНЦУ описывается кинетическим уравнением:

где $\alpha $ и $0 < \gamma < 1$ – параметры модели, определяющие скорость процесса развития усталостных повреждений. Выбор знаменателя в этом двухпараметрическом уравнении задает бесконечно большую скорость роста зоны полного разрушения при $\psi \to 1$ и обусловлен известными экспериментальными данными по кинетическим кривым роста усталостных трещин, имеющим вертикальную асимптоту и отражающим факт их взрывного, неконтролируемого роста на последней стадии макроразрушения.

Уравнение для повреждаемости сходного типа рассматривалось ранее [21, 22], его многочисленные параметры и коэффициенты определялись косвенно по результатам одноосных усталостных испытаний. Определим коэффициент ${{B}_{{LH}}}$ по процедуре, явно ассоциированной с выбранным критерием многоосного усталостного разрушения, которая выглядит следующим образом. Интегрируя уравнение для поврежденности при однородном напряженном состоянии

находим число циклов до полного разрушения ${{N}_{{LH}}}$ при $\psi = 1$

Приравнивая значения ${{N}_{{LH}}}$ из критерия разрушения (1.2) и из решения уравнения для повреждаемости (1.3) , получим выражение для коэффициента ${{B}_{{LH}}}$:

где величина ${{\sigma }_{{LH}}}$ определяется выбранным механизмом усталостного разрушения и соответствующим многоосным критерием (1.1).

2. Кинетическое уравнение для повреждаемости в режиме СВМУ.

2.1. Критерий многоосного усталостного разрушения в режиме СВМУ, соответствующий правой ветви бимодальной усталостной кривой на рис. 2 (обобщенный stress-based SWT) имеет вид:

где из условия подобия опорных точек для левой и правой ветвей бимодальной усталостной кривой [24] можно получить формулу ${{\sigma }_{V}} = {{10}^{{8{{\beta }_{{VH}}}}}}({{\sigma }_{u}} - {{\tilde {\sigma }}_{u}})$.

Из критерия усталостного разрушения получаем число циклов до разрушения при однородном напряженном состоянии:

Здесь ${{\tilde {\sigma }}_{u}}$ – предел усталости материала при реверсивном цикле для режима СВМУ, ${{\beta }_{{VH}}}$ – степенной показатель левой ветви бимодальной усталостной кривой.

2.2. Для режима СВМУ подобно п. 1.4 можно определить коэффициент в эволюционном уравнении для повреждаемости:

Аналогично уравнению (1.3) имеем:

Приравнивая значения ${{N}_{{VH}}}$ из критерия разрушения (2.2) и из решения уравнения для повреждаемости (2.3), получим выражение для коэффициента ${{B}_{{VH}}}$ в уравнении для повреждаемости:

где величина ${{\sigma }_{{VH}}}$ определяется выбранным механизмом усталостного разрушения и соответствующим ему многоосным критерием (2.1).

Поведение поврежденного материала описывается уравнениями деградации его модулей упругости $\lambda (\psi ) = {{\lambda }_{0}}(1 - \psi )$, $\mu (\psi ) = {{\mu }_{0}}(1 - \psi )$. Материал становится полностью разрушенным, его модули упругости обнуляются при достижении функции повреждаемости $\psi (N)$ значения 1.

3. Условие переключения режимов накопления усталостной поврежденности. Точка перехода с левой ветви усталостной кривой на правую ветвь, при которой происходит смена механизма усталостного разрушения, находится чуть выше предела усталости ${{\sigma }_{u}}$ (рис. 1) и определяется значением ${{\sigma }_{*}} = {{\sigma }_{u}} + \Delta \sigma $. Для обеспечения непрерывного сопряжения левой и правой ветвей усталостной кривой необходимо выполнить условие ${{N}_{{LH}}} = {{N}_{{VH}}}$, которое эквивалентно уравнению для величины $\Delta \sigma $:

или

С учетом фактической малости поправочного члена в квадратных скобках, можно установить значение поправки $\Delta \sigma $ приближенной формулой:

Соответствующее приближенное значение ${{N}_{*}} = {{N}_{{LH}}}({{\sigma }_{*}})$ определяется величиной ${{N}_{*}} = {{10}^{3}}{{\left[ {({{\sigma }_{B}} - {{\sigma }_{u}}){\text{/}}\Delta \sigma } \right]}^{{1/{{\beta }_{{LH}}}}}} \approx {{10}^{8}}$, как и следовало ожидать.

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

Для режима МЦУ–МНЦУ при ${{\sigma }_{u}} + \Delta {{\sigma }_{u}} < {{\sigma }_{{LH}}} < {{\sigma }_{B}}$, $\Delta \sigma = {{10}^{{ - 5{{\beta }_{{LH}}}}}}({{\sigma }_{B}} - {{\sigma }_{u}})$

Для режима СВМУ при ${{\tilde {\sigma }}_{u}} < {{\sigma }_{{VH}}} \leqslant {{\sigma }_{u}} + \Delta {{\sigma }_{u}}$

При ${{\sigma }_{{VH}}} \leqslant {{\tilde {\sigma }}_{u}}$ усталостного разрушения не происходит, при ${{\sigma }_{{LH}}} \geqslant {{\sigma }_{B}}$ оно наступает мгновенно.

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

Для расчета цикла нагружения деформируемого образца с дефектами был использован программный пакет Астра, использующий высокоэкономичный безматричный вариант МКЭ [28] и дополненный кодом для расчета уравнения усталостной повреждаемости и изменения модулей упругости.

Ранее этот вариант МКЭ применялся для развития методов сквозного счета процессов динамического разрушения [29, 30]. С его использованием были построены [31] алгоритмы решения контактных задач. Он применялся [32] для решения задач прессования и спекания образцов при квазистатических или нестационарных (нециклических) нагрузках.

Для интегрирования уравнения $d\psi {\text{/}}dN = B{{\psi }^{\gamma }}{\text{/}}(1 - {{\psi }^{\alpha }})$, $B = {{B}_{{LH}}}$, ${{B}_{{VH}}}$ применялась аппроксимация функции повреждаемости в k-узле расчетной сетки при заданных дискретных значениях $\psi _{k}^{n}$ в моменты ${{N}^{n}}$ и искомых $\psi _{k}^{{n + 1}}$ в моменты ${{N}^{{n + 1}}}$.

Отметим, что вблизи состояния полного разрушения при $\psi \to 1$ знаменатель кинетического уравнения становится сколь угодно малым, уравнение становится “жестким” и обычные явные методы его численного решения становятся непригодными. При численном решении по неявным схемам возникает нелинейная алгебраическая система, решение которой при большом числе шагов затратно по времени расчета и достижению приемлемой точности.

Поэтому был выбран следующий метод расчета повреждаемости. Было замечено, что для значения $\alpha = 1 - \gamma $ путем аналитического интегрирования кинетического уравнения для повреждаемости на приращении числа циклов можно получить явное выражение для $\psi _{k}^{{n + 1}}(\psi _{k}^{n},\Delta {{N}^{n}})$. Кинетическое уравнение при этом становится однопараметрическим, но сохраняет желаемые особенности поведения скорости роста повреждаемости (заданный параметром $\gamma $ степенной рост в окрестности нуля, и бесконечный рост при $\psi \to 1$).

Выполним это интегрирование:

Обозначив ${{(\psi _{k}^{{n + 1}})}^{{1 - \gamma }}} = x$, $q = 2(1 - \gamma )B\Delta {{N}^{n}} + {{(\psi _{k}^{n})}^{{1 - \gamma }}} - 2{{(\psi _{k}^{n})}^{{2(1 - \gamma )}}}$, $\Delta {{N}^{n}} = {{N}^{{n + 1}}} - {{N}^{n}}$, получим уравнение ${{x}^{2}} - 2x + q = 0$ и его допустимый корень $x = 1 - \sqrt {1 - q} < 1$. Из него следует точная формула для повреждаемости на приращении числа циклов $\Delta {{N}^{n}}$:

Значение приращения $\Delta {{N}^{n}}$ определяется следующим образом.

На основе данных расчета напряженного состояния для каждого узла рассчитывается коэффициент В. Для определения приращения числа циклов выбираем точку с наиболее быстро растущей функцией поврежденности, то есть ту, в которой раньше других $\psi _{k}^{n}$ возрастет на заданную величину $\Delta {{\psi }_{0}}$. Выборка проводится по точкам, где $0 \leqslant \psi _{k}^{n} \leqslant 1$. Величину шага определяем следующим образом:

Величина шага ограничена заданным максимумом $\Delta {{N}_{{\max }}}$. В расчетах используются значения $\Delta {{\psi }_{0}} = 0.1$, $\Delta {{N}_{{\max }}} = {{10}^{5}}$.

Для каждого узла, исходя из его текущего уровня повреждаемости и эквивалентного напряжения, находится новый уровень повреждаемости с учетов вычисленного приращения $\Delta {{N}^{n}}$.

Деградация поврежденного материала при численных расчетах определялась зависимостью $E_{k}^{{n + 1}} = {{E}_{0}}(1 - \psi _{k}^{{n + 1}})(H({{\psi }_{0}} - \psi _{k}^{{n + 1}})$ + 0.001), где $E_{k}^{{n + 1}}$ – значение модуля Юнга на n + 1-м шаге по циклам нагружения при фиксированном значении коэффициента Пуассона, H – функция Хэвисайда, а ${{\psi }_{0}}\sim 1$ – критический параметр повреждаемости, по достижении которого материал переходит в полностью разрушенное состояние с околонулевыми модулями упругости. При расчетах выбиралось ${{\psi }_{0}} = 0.9$. Малое значение остаточного модуля Юнга – 0.001E для полностью разрушенного материала играет роль регуляризатора и вводится для того, чтобы сохранялась возможность вести сквозной счет для любого состояния поврежденности.

Расчет заканчивается при выходе границ области полностью разрушенного материала на поверхность образца (макроразрушение) или прекращении эволюции этой области.

С помощью мультирежимной модели и разработанной численной процедуры решены задачи о развитии зон усталостной поврежденности в окрестности типичных дефектов структуры материала (микропоры, различным образом ориентированные микротрещины).

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

5. Примеры расчетов. Были проведены расчеты циклического нагружения пластин с дефектами для коэффициента асимметрии R = –1 (реверсивный цикл) c зарождением зон поврежденности, возникновения и развития трещиноподобных зон полного разрушения частиц материала (“квазитрещин”) вплоть до их выхода на боковую грань пластины (макроразрушения).

Нормальная нагрузка прикладывается с торцов, боковые грани пластины свободны от напряжений. Рассматривается плоское напряженное состояние. Упругие, прочностные и усталостные свойства материала пластины (титановый сплав) были взяты из [21]: ${{\sigma }_{B}} = {\text{116}}0$ МПа, ${{\sigma }_{u}} = 337$ МПа, ${{\tilde {\sigma }}_{u}} = 250$ МПа, ${{\beta }_{{LH}}} = 0.31$, ${{\beta }_{{VH}}} = 0.27$. Модули упругости неповрежденного сплава равны ${{\lambda }_{0}} = 77$ ГПа, ${{\mu }_{0}} = 44$ ГПа.

Степенной показатель кинетического уравнения $\gamma $ должен быть определен в процессе согласования результатов вычислительных и усталостных экспериментов. Для численных расчетов в уравнении для повреждаемости было принято среднее значение $\gamma $ = 0.5, $\alpha = 1 - \gamma $.

Размеры пластины 40 × 40 мм. Горизонтальная y = 0 и вертикальная x = 0 центральные оси пластины являются осями симметрии, поэтому расчет проводится для области x > 0, y > 0. Рассмотрены три вида дефектов в центре пластины: круговое отверстие радиусом r = 1 мм; узкое эллиптическое отверстие с горизонтальной осью a = 1 мм, вертикальной осью b = 0.2a; узкое эллиптическое отверстие с вертикальной осью b = 1 мм, горизонтальной осью a = 0.2b. На рис. 2 показаны результаты расчета для кругового отверстия. На этом и последующих рисунках представлены изолинии максимального главного напряжения вокруг “квазитрещины” в тот момент, когда она прошла примерно половину расстояния до боковой грани. В подписях к рисункам указаны значения количества циклов: ${{N}_{i}}$ – до начала разрушения, ${{N}_{f}}$ – до макроразрушения, ${{N}_{k}}$ – для положения, представленного на графике. Также указан максимальный уровень растягивающих напряжений Р в цикле приложенных нагрузок и соответствующий им режим усталостного разрушения.

На рис. 3 показаны результаты расчета для горизонтально ориентированного узкого эллиптического отверстия; на рис. 4 – для вертикально ориентированного узкого эллиптического отверстия.

Отметим, что количество циклов ${{N}_{k}}$, соответствующее уровню продвижения трещиноподобных зон разрушения частиц материала (“квазитрещин”) примерно до середины образца (рис. 2–4), очень мало отличается от количества циклов ${{N}_{f}}$, что указывает на дальнейший финальный, быстропротекающий процесс макроразрушения.

Из рисунков видно, что близкий уровень циклических нагрузок, приложенных к образцу с дефектами разного типа и ориентации, может приводить к усталостному макроразрушению образца в совершенно различных режимах. Так, например, для вертикально ориентированного эллиптического дефекта при уровне P = 264 МПа имеем режим, переходный к СВМУ (рис. 4а) с долговечностью ${{N}_{f}} = 4.709 \times {{10}^{7}}$, а для горизонтально ориентированного эллиптического дефекта при существенно более низком P = 180 МПа макроразрушение происходит по механизму МНЦУ (рис. 3б) с долговечностью ${{N}_{f}} = 1.057 \times {{10}^{6}}$.

Также можно отметить, что при проведении расчетов по мультирежимной модели неоднократно проявлялся и моделировался эффект, наблюдаемый при фрактографических исследованиях элементов конструкций, разрушенных при эксплуатации [33]. При относительно небольшом уровне приложенных нагрузок, за счет внутренних дефектов зарождение и начальное развитие зон усталостного разрушения происходит по механизму СВМУ, а затем при концентрации и росте уровня напряжений в окрестности фронта распространяющейся “квазитрещины” происходит переход на другую ветвь (МЦУ–МНЦУ) усталостной кривой. При этом резко меняется скорость развития процесса усталостного разрушения.

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

Заключение. Предложена мультирежимная кинетическая модель развития повреждаемости при циклическом нагружении для описания развития процесса усталостного разрушения. Для определения коэффициентов кинетического уравнения повреждаемости использован известный критерий многоосного усталостного разрушения SWT, в котором заложен механизм, связанный с развитием микротрещин нормального отрыва. На этой основе предложена процедура вычисления коэффициентов кинетического уравнения для различных режимов усталостного разрушения МЦУ–МНЦУ, СВМУ.

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

Исследование выполнено за счет гранта Российского научного фонда (проект № 19-19-00705).