Fatigue behavior under high frequency loading of materials produced by selective laser melting

Full Text

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

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

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

При определенном увеличении мощности луча или его замедлении и при достаточном перекрытии соседних по горизонтали и вертикали треков возникают периодические системы зон переплавов. Такая структура, несмотря на отсутствие пор, также является в какой-то мере дефектной, поскольку при дополнительном переплавлении отвержденных зон меняется строение кристаллической решетки сплава. Исходная альфа-фаза титанового сплава при переплавлении переходит в бета-фазу, а при дополнительном переплавлении в альфа+бета [2]. Эти состояния имеют отличающиеся механические характеристики, по измерениям микротвёрдости – до 20%. Упругие характеристики также могут отличаться на 5–10% [2, 3]. Таким образом, материал, полученный с использованием СЛП технологии, в самом благоприятном случае имеет слабо неоднородную периодическую структуру, определяемую расположением зон многократных переплавов.

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

Следует отметить, что диапазон параметров мощности-скорости лазерного луча для создания относительно благоприятной структуры СЛП материала является достаточно узким и его определение расчетным образом представляет собой отдельную серьезную задачу.

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

Статические прочностные испытания образцов из СЛП сплавов показывают, что их упругие модули и статический предел прочности мало отличаются от макроскопических механических характеристик образцов, полученных традиционными технологиями [5]. Однако высокочастотные (20 кГц) циклические испытания на пьезоэлектрической установке Lazur [6, 7] показали, что в режиме сверхмногоциклового усталостного разрушения (СВМУ) [8] долговечность СЛП образцов может падать на один-два порядка для материалов с сильно развитой системой дефектов типа недоплавленных пор. Но даже для образцов со слабо неоднородной структурой, возникающей за счет периодической системы переплавленных зон, долговечность при таких испытаниях уменьшается в 2–3 раза [6, 7]. Этот эффект можно объяснить более быстрым зарождением начальных повреждений на внутренних структурных границах за счет повышенной концентрации напряжений, а также повышенной чувствительностью количества циклов до разрушения к уровню амплитуды напряжений в высокочастотном цикле в режиме СВМУ за счет более пологого ниспадания правой ветви бимодальной усталостной S-N кривой [8, 9].

Для исследования усталостной прочности СЛП образцов в режиме СВМУ при длительных высокочастотных нагрузках будет использована предложенная ранее мультирежимная модель циклической повреждаемости [9, 10]. Разработанные на ее основе вычислительные алгоритмы позволяют рассчитать долговечность образцов для различных систем дефектов микроструктуры, соответствующих заданным характеристикам подвижного лазерного луча, а также определить диапазон параметров процесса СЛП, при котором будут достигаться наилучшие показатели усталостной прочности и долговечности.

2. Моделирование процесса селективного лазерного плавления. Графическая схема постановки задачи о подвижном тепловом потоке от лазерного луча при его воздействии на слой порошка с образованием ванны расплава представлена на рис. 1.

 

Рис. 1. Графическая схема задачи о подвижном тепловом потоке от лазерного луча 1 при его воздействии на слой порошка 2 с образованием ванны расплава 6. Цифрой 3 обозначена подвижная фазовая граница, 4 – подложка из материала в твердой фазе, 5 – остывший след лазера в виде материала в твердой фазе.

 

Для математической постановки задачи селективного лазерного плавления будем использовать энтальпийную формулировку многофазной нелинейной задачи теплопроводности, не требующую явного выделения неизвестных межфазных границ [11, 12], и в последние годы используемую различными авторами [13–16].

Эта формулировка выглядит следующим образом.

Уравнение для энтальпии:

$\frac{\partial E}{\partial t}=\frac{\partial }{\partial x}\left(k\left(E\right)\frac{\partial T\left(E\right)}{\partial x}\right)+\frac{\partial }{\partial y}\left(k\left(E\right)\frac{\partial T\left(E\right)}{\partial y}\right)+\frac{\partial }{\partial z}\left(k\left(E\right)\frac{\partial T\left(E\right)}{\partial z}\right).$ (2.1)

Граничные условия на гранях параллелепипеда имеют следующий вид.

Верхняя грань z = 0:

$\begin{array}{c}{\left.k\left(E\right)\frac{\partial T}{\partial z}\right|}_{z=0}=-b\left(T_{z=0}-T_{ext}\right)+q(x-vt,y)/\left(p{R}^2\right),\\ q(x-vt,y)=q_{0}H\left(R-r\right),\end{array}$ (2.2)

где R – радиус пятна лазерного луча, q₀ – полная мощность лазерного луча, $r\left(x-vt,y\right)=\sqrt{{(x-vt)}^2+y^2}$ – расстояние до центра лазерного пятна, который движется со скоростью v вдоль оси x, H(x) – функция Хэвисайда.

На остальных гранях:

$\begin{array}{cc}x=\mathrm{0,}x=a& k\left(E\right){\left.\frac{\partial T}{\partial x}\right|}_{x=\mathrm{0,}a}\text{}=-b\left(T_{x=\mathrm{0,}a}-T_{ext}\right)\\ y=-b/\mathrm{2,}y=b/2& k\left(E\right){\left.\frac{\partial T}{\partial y}\right|}_{y=-b/\mathrm{2,}b/2}\text{}=-b\left(T_{y=-b/\mathrm{2,}b/2}-T_{ext}\right)\\ z=-(h+d)& k\left(E\right){\left.\frac{\partial T}{\partial z}\right|}_{z=-(h+d)}\text{}=-b\left(T_{z=-(h+d)}-T_{ext}\right)\end{array}$ (2.3)

Начальные условия:

${\left.E(x,y,z,t)\right|}_{t=0}=E\left(T_0\right)={\rho }_0(x,y,z)c_0(x,y,z)T_0,$

где p₀(z), c₀(z), k₀(z) – начальные характеристики плотности, теплоемкости и теплопроводности порошковой и твердой фазы, определяемые вертикальной координатой точки расположения, T_ext – температура внешней среды, b – коэффициент теплопередачи «металл – внешняя среда», например: при $z<-h,{\rho }_0\left(z\right)={\rho }_S,$ при $z>-h,{\rho }_0\left(z\right)={\rho }_P.$ Здесь и далее индексами S и P помечены значения параметров порошковой и твердой фаз.

Принимается следующая зависимость температуры от энтальпии, учитывающая возможные фазовые переходы:

$\mathit{T}\mathbf{=}\{\begin{array}{ll}\frac{E}{r_0\left(z\right)c_0\left(z\right)},& E<E_{m-},\\ T_m,& E_{m-}\le E<E_{m+},\\ T_m+\frac{(E-E_{m+})}{r_L{c}_L},& E_{m+}\le E<E_{v-},\\ T_v,& E_{v-}\le E<E_{v+},\\ T_v+\frac{(E-E_{v+})}{r_V{c}_V},& E_{v+}\le E,\end{array}$ (2.4)

$E_{m-}={\rho }_0\left(z\right)c_0\left(z\right)T_m,$ $E_{m+}=E_{m-}+{\rho }_0\left(z\right)l,$ $E_{v-}={\rho }_L{c}_L(T_v-T_m)+E_{m+},$ E_v+ = E_v–+ +p_VL, где p_L, c_L и k_L – характеристики (плотность, коэффициенты теплоемкости и теплопроводности) жидкой фазы, T_m – температура плавления, $\lambda$ – удельная теплота плавления, p_L, c_L, k_L – характеристики (плотность, коэффициенты теплоемкости и теплопроводности) газообразной фазы, T_v – температура парообразования, $\Lambda$– удельная теплота парообразования.

Зависимость теплопроводности от энтальпии с учетом фазовых переходов имеет вид:

$k\left(E\right)=\left\{\begin{array}{ll}{k}_0\left(z\right),& E<E_{m-},\\ k_0\left(z\right)+\left[k_L-k_0\left(z\right)\right]\frac{(E-E_{m-})}{(E_{m+}-E_{m-})},& E_{m-}\le E<E_{m+},\\ k_L,& E_{m+}\le E<E_{v-},\\ k_L+\left[k_V-k_L\right]\frac{(E-E_{v-})}{(E_{v+}-E_{v-})},& E_{v-}\le E<E_{v+},\\ k_V& E_{v+}\le E,\end{array}\right.$ (2.5)

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

3. Численное моделирование процесса СЛП.

3.1. Численная схема решения многофазной задачи о формировании ванны расплава и отвержденного трека. Для аппроксимации энтальпийного уравнения (2.1) с заданными граничными и начальными условиями используется явная разностная схема в варианте метода конечных объемов [17]. Аппроксимация во внутренних точках основана на интегрировании уравнения (2.1) по объему ячейки с центральной точкой определения неизвестных сеточных функций $E_{ijk}^{n+1},$ $T_{ijk}^{n+1}$ на слое n + 1 по времени и имеет вид:

$\begin{array}{c}\underset{V}{\iiint }\frac{\partial E}{\partial t}dv=\underset{S}{∯}\left(Q_x{n}_x+Q_y{n}_y+Q_z{n}_z\right)ds,\\ Q_x=k\left(E\right)\frac{\partial T\left(E\right)}{\partial x},Q_y=k\left(E\right)\frac{\partial T\left(E\right)}{\partial y},Q_z=k\left(E\right)\frac{\partial T\left(E\right)}{\partial z},\\ \frac{E_{ijk}^{n+1}-E_{ijk}^n}{\Delta t}\Delta v_{ijk}=\sum _l{\left(Q_x{n}_x+Q_y{n}_y+Q_z{n}_z\right)}_{ijk}^l\Delta s_{ijk}^l,\end{array}$

где l – нумерация граней ячейки, x_i, y_j, z_k – координаты центров ячеек, на которые разбивается область решения, i = 1..IM, j = 1..JM, k = 1..KM.

В рассматриваемом случае внутренних ячеек-параллелепипедов с гранями, параллельными осям декартовой системы координат, тепловые потоки через грани ячейки равны:

$\begin{array}{c}{\left(Q_x{n}_x\right)}_{ijk}={\left(Q_x^{+}\right)}_{ijk}-{\left(Q_x^{-}\right)}_{ijk},\\ {\left(Q_x^{+}\right)}_{ijk}=\frac{(k_{i+\mathrm{1,}jk}^n+k_{ijk}^n)}{2}\frac{(T_{i+\mathrm{1,}jk}^n-T_{ijk}^n)}{(x_{i+\mathrm{1,}jk}-x_{ijk})},\\ {\left(Q_x^{-}\right)}_{ijk}=\frac{(k_{i-\mathrm{1,}jk}^n+k_{ijk}^n)}{2}\frac{(T_{i,jk}^n-T_{i-\mathrm{1,}jk}^n)}{(x_{ijk}-x_{i-\mathrm{1,}jk})},\end{array}$ (3.1)

и аналогичные выражения для направлений по осям x и y.

Формулы для значений энтальпии во внутренних ячейках в итоге имеют вид:

$\begin{array}{c}{E}_{ijk}^{n+1}=E_{ijk}^n+\sum _l{\left(Q_x{n}_x+Q_y{n}_y+Q_z{n}_z\right)}_{ijk}^l\Delta s_{ijk}^l\Delta t/\Delta v_{ijk},\\ i=\mathrm{2..}IM-\mathrm{1,}j=\mathrm{2..}JM-\mathrm{1,}k=\mathrm{2..}KM-1.\end{array}$ (3.2)

При этом учитываются принятые зависимости (2.4 – 2.5) для  $T_{ijk}^n\left(E_{ijk}^n\right),$ $k_{ijk}^n\left(E_{ijk}^n\right).$

Аппроксимация граничных ячеек на примере верхней грани по форме совпадает с (3.2) с учетом значения k = KM для граничной ячейки на верхней грани, но значение потока вычисляется с учетом его граничного значения:

$Q_z^{+}=-\beta \left(T_{ijk}^b-T_{ext}\right)+q(x-vt,y).$

Шаг по времени выбирается из условия устойчивости для выбранной явной аппроксимации трехмерного нестационарного параболического уравнения и имеет стандартный вид:

$\Delta t=\underset{ijk}{\min }\left(h^2\right)/\underset{ijk}{\max }(k/(\rho c\left)\right)/8.$

3.2. Результаты расчетов для единичного трека и сравнение с экспериментом. Были проведены расчеты единичного затвердевшего трека и ванны расплава в порошковом слое на подложке из титанового сплава Ti-6Al-4V для подвижного лазерного луча с мощностью $q_0=100 Вт,$ скоростью движения $\nu =800 мм/с,$ радиусом пятна $R=40 мкм.$ Расчеты проводились на сетке со сгущением в окрестности ванны расплава и минимальным пространственным шагом 2.5 мкм. Физические свойства трёх состояний, представленных в задаче и использованных при расчете – твёрдого титанового сплава, порошка из него и расплава – представлены в таблице 1.

 

Таблица 1. Характеристики Ti-6Al-4V и AlSi10Mg

Характеристика			Величина
Название	Обозначение	Размерность	Ti-6Al-4V	AlSi10Mg
Плотность твёрдой фазы	${р}_S$	кг/м3	4400	2670
Насыпная плотность порошка	${р}_P$	кг/м3	3750	2100
Плотность жидкой фазы	${р}_L$	кг/м3	3900	2460
Температура плавления	Tm	К	1920	870
Температура испарения	TV	К	3530	2600
Удельная теплота плавления	$\lambda$	кДж/кг	286	440
Удельная теплота парообразования	L	МДж/кг	4.0	8.0
Теплопроводность твёрдой фазы, Т = 300 К	kS	Вт/(м · К)	3.9	148.3
Теплопроводность твёрдой фазы, T = Tm	kS	Вт/(м · К)	28.1	136.5
Теплопроводность жидкой фазы, T = Tm	kL	Вт/(м · К)	28.6	77.4
Теплопроводность жидкой фазы, T = TV	kL	Вт/(м · К)	57.7	118.0
Теплопроводность газа, T = TV	kV	Вт/(м · К)	60.0	120.0
Теплопроводность порошка	kP	–	0.95	0.90
Теплоёмкость твёрдой фазы	CS	Дж/кг/К	760	990
Теплоёмкость порошка	CP	Дж/кг/К	760	990
Теплоёмкость жидкой фазы	CL	Дж/кг/К	1150	1470
Коэффициент теплопередачи металл-воздух	$\beta$	Вт/м2/К	10–100	10–100

 

На рис. 2 представлены результаты расчета формы ванны расплава и единичного затвердевшего трека. Различным цветом выделены фазовые состояния материала (синим – порошок, зелёным – расплав, красным – твёрдый сплав). Хорошо видно, что при данных параметрах процесса происходит желаемое сплавление трека с подложкой, но с достаточно малой зоной сцепления, составляющей примерно 20% от общей глубины трека.

 

Рис. 2. Одиночный трек СЛП в титановом сплаве, q₀ = 100 Вт, v = 800 мм/с, вид сверху – a), вид на сечение трека и его сплавление с подложкой – b), где w – ширина трека, d – глубина.

 

На основе серии расчетов были получены данные о ширине и глубине сечений сформировавшихся треков и построены их зависимости от скорости сканирования лазерного луча в диапазоне 500–1200 мм/c и для различных значений мощности из достаточно широкого диапазона 50-200 Вт (рис. 3, закрашенные точки). Эти зависимости сравнивались с результатами экспериментов из [1] (рис. 3, полые точки). Для получения этих результатов в зависимостях (2.4)–(2.5) для T(E) и k(E) был учтен фазовый переход жидкость-пар при температуре парообразования T_v. При этом учитывались дополнительные энергозатраты на удельную теплоту парообразования $\Lambda$. В зонах превышения температуры T_v, которые возникали при высоких значениях мощности лазерного луча, применялись дополнительные коэффициенты теплопроводности k_V и теплоемкости при значениях $\Lambda =4$МДж, k_V = 1200 Вт/(м · °С), c_V = 1150 Дж/кг/°С. Без учета этого фактора расхождение численных и экспериментальных результатов было существенно выше.

 

Рис. 3. Результаты расчётов ширины w – а) и глубины d – b) трека при различных скоростях движения лазерного луча v (мм/с) с учетом испарения с поверхности ванны расплава. Экспериментальные данные «1» из работы [1] сравнены с результатами расчета «2».

 

3.3. Результаты для линейной и многослойной системы треков.

3.3.1. Результаты расчетов для алюминиевого сплава. Микроскопическая структура напечатанного материала видна на снимках поверхности (выделенный трек, рис. 4, а) и среза сечения (система перекрывающихся треков, рис. 4, b), полученных с помощью электронного и оптического микроскопов при различных сочетаниях мощности и скорости движения лазерного луча.

 

Рис. 4. Наблюдаемая в электронный микроскоп структура напечатанного алюминиевого сплава, а) – вид сверху, q₀ = 370 Вт, v = 650 мм/с, электронный микроскоп, b) – сечение приповерхностного слоя, q₀ = 370 Вт, v = 450 мм/с, оптический микроскоп.

 

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

 

Рис. 5. Зоны переплавов (синий цвет – 0 переходов, жёлтый – 3 перехода, оранжевый – 4 перехода) для трехслойной печати – а), четырехслойной печати – b).

 

Для трехслойной (рис. 5, a) и четырехслойной (рис. 5, b) печати разным цветом показаны области с многократным переходом из жидкого в твердое состояние (переплавы, синий цвет – 0 переходов, жёлтый – 3 перехода, оранжевый – 4 перехода).

Было проведено прямое численное моделирование процесса создания подложки в алюминиевом порошке и сравнение с реально напечатанной подложкой (рис. 6). Для создания подложки использовалось селективное лазерное плавление со следующими параметрами процесса: диаметр пятна 80 мкм, мощность 320 Вт, скорость сканирования 1500 мм/с, шаг луча 105 мкм. Можно сказать, что геометрическая структура перекрывающихся треков воспроизводится достаточно хорошо.

 

Рис. 6. Геометрическая структура перекрывающихся треков. Верхняя часть – расчет, нижняя часть – эксперимент, оптический микроскоп.

 

3.3.2. Результаты расчетов для титанового сплава. Аналогичные расчеты для определения структуры дефектов различного типа были проведены для титанового сплава. На рис. 7, а показана смешанная структура дефектов, возникающая при определенных параметрах процесса СЛП, q₀ = 50 Вт, v = 1000 мм/с. В начале процесса видны недоплавы (дефекты в виде пор, синий цвет), которые постепенно уходят при многократном многослойном прохождении лазерного луча и трансформируются в периодическую структуру неоднородностей (переплавов, желтый и оранжевый цвет). На рис. 7, b при параметрах q₀ = 50 Вт, v = 500 мм/с, с уменьшением скорости движения луча, видны множественные переплавы (желтый и оранжевый цвет), недоплавы (поры) отсутствуют.

 

Рис. 7. Результаты расчета структуры дефектов для титанового сплава, a) – q₀ = 50 Вт, v = 1000 мм/с, недоплавы, b) – q₀ = 50 Вт, v = 500 мм/с, множественные переплавы, недоплавы отсутствуют.

 

Таким образом, послойное переплавление перекрывающихся треков создает регулярную, квазипериодическую неоднородную структуру с характерными размерами, по порядку величины совпадающую с поперечными размерами треков, т.е. 100–500 мкм в зависимости от параметров процесса СЛП. Размеры дефектов при этом составляют 30–40 мкм.

Эти данные использованы для дальнейшей оценки наблюдаемого эффекта падения усталостной прочности СЛП образцов [6, 7] с помощью модели циклической повреждаемости при длительных высокочастотных испытаниях [9, 10].

4. Модель циклической повреждаемости образцов при высокочастотном нагружении.

4.1. Формулировка мультирежимной двухкритериальной модели. Развитие представлений об усталостном разрушении материалов [18] связано с экспериментальными исследованиями процессов циклического деформирования с характерной долговечностью в диапазоне 10³ – 10⁵ циклов (режим малоцикловой усталости – МЦУ), и циклическом нагружении, приводящем к разрушению в диапазоне 10⁵ – 10⁷ циклов (многоцикловая усталость – МНЦУ). Ранее считалось, что эти два диапазона составляют полную кривую усталости с выходом на «предел усталости», ниже которого усталостного разрушения не происходит. С развитием средств управления нагружающими машинами стало возможно проводить высокочастотные усталостные испытания при значительно больших долговечностях (порядка 10⁸–10¹⁰ циклов) [19]. Оказалось, что при таких наработках материалы могут разрушаться при уровнях амплитуд напряжений существенно ниже классического «предела усталости» (режим сверхмногоцикловой усталости – СВМУ) [20]. Зарождение трещины в этом режиме связано с микроскопическими особенностями строения материала или наличием внутренних дефектов и происходит под поверхностью образца. Результаты испытаний в МЦУ и МНЦУ режимах усталостного разрушения в основном отражают макромасштабные свойства сплавов, тогда как высокочастотные испытания в СВМУ режиме отражают микромасштабное усталостное поведение. Поскольку при процедуре 3Д-печати по технологии селективного лазерного плавления в образце создается довольно сложная микроструктура с различными особенностями, СВМУ испытания дают хорошую возможность исследования микромасштабного усталостного поведения СЛП материалов.

Сформулируем связанную задачу о высокочастотном нагружении образца, в котором могут развиваться процессы циклической повреждаемости, описываемые функцией повреждаемости y [10], для случая стационарных гармонических колебаний: w(x, t) = u(x)e ^iwt.

Уравнения движения и закона Гука для стационарных колебаний неоднородного образца имеют вид:

$\begin{array}{c}\nabla \times s+\rho {\omega }^{2}u=0,\\ s=\lambda \left(\psi \right)(e:I)I+2\mu \left(\psi \right)e,\\ \epsilon =\left(\nabla \otimes u+{(\nabla \otimes u)}^T\right)/2.\end{array}$ (4.1)

Система уравнений в каждой точке границы дополняется либо кинематическими граничными условиями:

$\begin{array}{c}{\left.u\cdot {\tau }_{\alpha }\right|}_{x\in \partial V_{u\alpha }}=u_{\alpha }^{*}(,t),\\ {\left.u\cdot n\right|}_{x\in \partial V_{un}}=u_n^{*}(,t)\end{array}$ (4.2)

либо динамическими граничными условиями:

$\begin{array}{l}{\left.:n{\otimes }_{\alpha }\right|}_{\in \partial V_{\tau \alpha }=\partial V/\partial V_{u\alpha }}=p_{\alpha }^{*}(,t),\\ {\left.:n\otimes n\right|}_{\in \partial V_n=\partial V/\partial V_{un}}=-p_n^{*}(,t),\end{array}$

где n и t_a(a = 1, 2) – орты нормали и касательных к границе. Заданные на границе функции отмечены звездочками. Модули упругости Ламе уменьшаются в зависимости от накопленной распределенной повреждаемости.

Изменение функции повреждаемости $\psi$ с ростом циклов нагружения N описывается кинетическим уравнением, предложенным в [9, 10]:

$d\psi /dN=B(\sigma ,\Delta \sigma ){\psi }^{\gamma }\text{}/\text{}(1-{\psi }^{1-\gamma }),$ (4.3)

где $\beta \left(\sigma ,∆\sigma \right)$ – коэффициент, зависящий от напряжённо-деформированного состояния в цикле нагружения, $∆\sigma$ – размах циклической нагрузки, $\gamma$ – параметр, определяющий скорость накопления повреждаемости.

Выражения для коэффициентов имеют вид [9]:

$B=\left\{\begin{array}{ccc}{B}_L=\frac{{\left[⟨{\sigma }_{eq}-{\sigma }_u⟩/\left({\sigma }_B-{\sigma }_u\right)\right]}^{1/{\beta }_L}}{2\left(1-\gamma \right){10}^3}& при& {\sigma }_u+\Delta {\sigma }_u<{\sigma }_{eq}<{\sigma }_B,\\ B_V=\frac{{\left[⟨{\sigma }_{eq}-{\tilde{\sigma }}_u⟩/\left({\sigma }_u-{\tilde{\sigma }}_u\right)\right]}^{1/{\beta }_V}}{2\left(1-\gamma \right){10}^8}& при& {\tilde{\sigma }}_u<{\sigma }_{eq}\le {\sigma }_u+\Delta {\sigma }_u,\end{array}\right.$

где $\Delta {\sigma }_u={10}^{-5{\beta }_L}\left({\sigma }_B-{\sigma }_u\right)$ – ширина области бифуркации [9]. Выражения, стоящие в треугольных скобках, указывают на соответствующую ветвь усталостной кривой и определяются как ⟨ f ⟩ = f H( f ), где H( f ) – функция Хэвисайда. Величина g лежит в диапазоне 0 < $\gamma$ < 1.

Разрушение по механизмам отрыва связано с растягивающими компонентами напряжений и отражается в критерии SWT [21]:

${\sigma }_{eq}={\sigma }^n=\sqrt{⟨{\sigma }_{1\max }⟩\Delta {\sigma }_1/2},$

где $⟨{\sigma }_{1\max }⟩$ – величина максимального главного растягивающего напряжения и где $∆{\sigma }_1/2$ – амплитуда главного циклического напряжения.

Для описания развития микротрещин сдвигового типа используется критерий CSV [22], учитывающий особенности их формирования:

${\sigma }_{eq}={\sigma }^{\tau }=\sqrt{{(⟨\Delta {\sigma }_n⟩/2)}^2+3{(\Delta {\tau }_n/2)}^2},$

где $\Delta {\tau }_n/2$ – амплитуда максимального касательного напряжения, действующей на некоторой площадке с нормалью n, $(⟨\Delta {\sigma }_n/2⟩$ – амплитуда циклического напряжения в фазе растяжения $⟨\Delta {\sigma }_n⟩={\sigma }_{n_{\max }}H\left({\sigma }_{n_{\max }}\right)-{\sigma }_{n_{\min }}H\left({\sigma }_{n_{\min }}\right),$ действующего на этой же площадке.

Развитие повреждаемости в материальной частице приводит к эффективному уменьшению модулей упругости. В предложенном варианте модели – по кусочно-линейному закону следующего вида:

$\lambda ,\mu \left(\psi \right)=\left\{\begin{array}{ccc}{\lambda }_0,{\mu }_0\left(1-\kappa \psi \right)& при& \psi <{\psi }_{*},\\ 0& при& {\psi }_{*}\le \psi \le \mathrm{1,}\end{array}\right.$ (4.4)

где $\psi * < 1$ – критическое значение повреждаемости, при котором наступает состояние полного разрушения.

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

4.1.1. Нелокальная повреждаемость. При численном решении задач усталостного разрушения в континуальной постановке с использованием функций циклической повреждаемости включение выбранных критериев SWT или CSV определяется классическим пределом усталости s_u и СВМУ пределом усталости $\widetilde{{\sigma }_u}$. При этом, как отмечено во многих работах по континуальному разрушению [23–25], возникают проблемы обеспечения сходимости интегральных критических нагрузок и долговечности образцов до разрушения при измельчении сетки. Поэтому приведенная выше модель усталостного разрушения была дополнена уравнениями, которые при численном расчете расширенной системы обеспечивают сходимость расчетных критических нагрузок и числа циклов до разрушения при измельчении сетки. Такое расширение достигается с помощью введения так называемой нелокальной повреждаемости, основанной на осреднении функции повреждаемости в пространственной области с заданным малым линейным размером [24–26]. Его величина определяется линейным параметром структуры материала, который можно принять равным размеру зерна r₀.

Локальная модель повреждаемости регуляризуется в нелокальный вариант дополнительным осреднением. Для этого вводится нормированный оператор делокализации, который применяется либо к самой функции повреждаемости $\psi$, либо к некоторой связанной с ней функции $f = f \left(\psi \right)$. Этот оператор может иметь различные представления, в частности форму [26]:

$\begin{array}{cc}\tilde{f}\left(x\right)=\underset{r\pounds 1}{\int }f\left(x\right)\overline{L}(x,y)dy,& \overline{L}(x,y)=L_0(x,y)/\underset{r\pounds 1}{\int }{L}_0(x,z)dz,\\ L_0(x,y)={⟨1-r^m(x,y)⟩}^2,& r^2(x,y)={\left|x-y\right|}^2\text{}/\text{}{r_0}^2.\end{array}$ (4.5)

Делокализация сводится к процедуре осреднения внутри «колокола» с основанием 2r₀ и эффективной шириной, определяемой степенным параметром m.

Корректировка модулей упругости (4.4) в этом варианте модели происходит уже в зависимости от осредненной повреждаемости $\widetilde{\psi ,} т.е. \lambda =\left(\tilde{\psi }\right), \mu =\mu \left(\tilde{\psi }\right)$.

При численной реализации такой процедуры существенную роль играет расстояние между узлами ячеек сетки h. Если h > r₀, то значение повреждаемости не меняется. Если же h < r₀, то происходит осреднение повреждаемости по узлам, попадающим в «колокол». Такое осреднение приводит к стабилизации ширины квазитрещины на уровне порядка r₀ и обеспечивает сходимость долговечности даже при значительном измельчении сетки.

4.1.2. Численный метод решения связанной системы уравнений. Численный метод расчета зон повреждаемости состоит в расчете упругого напряженно-деформированного состояния образца по циклам нагружения, численного решения нелинейного уравнения (4.3) для повреждаемости $\psi$ и вычисления модулей упругости среды в зависимости (4.4) от $\psi$.

Разностная аппроксимация уравнения (4.3) выполняется путем интегрирования на интервале между двумя дискретными значениями циклов N ⁿ и N ^{n + 1}. Тогда для функции повреждаемости в каждом пространственном узле сетки с заданным шагом по циклам нагружения можно получить аналитическое решение:

${\psi }_k^{n+1}={\left(1-\sqrt{{\left(1-{\left({\psi }_k^n\right)}^{1-\gamma }\right)}^2-2\left(1-\gamma \right)B_k^n\Delta N^n}\right)}^{1/(1-\gamma )}\text{},$ (4.6)

где ${\psi }_k^{n+1}$ – значение функции повреждаемости в пространственном узле k и временном слое n + 1, и где $\Delta N^n=N^{n+1}-N^n\text{}.$
Шаг расчета по числу циклов определен в [9]:

$\Delta N^n=\underset{k}{\min }\Delta {\tilde{N}}_k^n/2,\Delta {\tilde{N}}_k^n={\left.\left(\frac{{\psi }^{1-\gamma }}{1-\gamma }-\frac{{\psi }^{2\left(1-\gamma \right)}}{2\left(1-\gamma \right)}\right)/B_k^n\right|}_{{\psi }_k^n}^1\text{}.$

Зависимость между значениями упругих характеристик и функции повреждаемости принята в виде (4.4), но с учетом описанного выше расширения модели (4.5) вычисляется в зависимости от нелокальной повреждаемости $\tilde{\psi } \left(х\right):$

$E_k^{n+1}=E_0\left(1-\kappa {\tilde{\psi }}_k^{n+1}\right)\left[H\left({\psi }_{*}-{\tilde{\psi }}_k^{n+1}\right)+0.001\right],$ (4.7)

где ${Е}_k^{n+1}$ – значение модуля упругости Юнга на новом шаге, E₀ – модуль Юнга неповреждённого материала, k – коэффициент деградации модуля упругости, определяется в ходе вычислительных экспериментов. В состоянии полного разрушения модуль Юнга принимается равным 0.001 от своего начального значения, что позволяет вести расчет сильно неоднородных упругих задач на каждом шаге по циклам нагружения для фиксированной сетки.

4.2. Примеры численных расчетов зарождения и развития повреждений для стандартных образцов. Схема расчета высокочастотных колебаний образцов с использованием модели циклической повреждаемости (4.1)–(4.5) и расчетных формул (4.6)–(4.7), была использована для моделирования их усталостного разрушения. Проверка сходимости численного решения при измельчении сетки с использованием процедуры делокализации была проведена в [10]. Расчеты проводились при следующих значениях параметров модели циклической повреждаемости, которые в среднем соответствуют характеристикам бимодальной усталостной кривой для стали при коэффициенте асимметрии R = -1 (симметричный цикл): E₀ = 207 ГПа,${\sigma }_B$ = 1600 МПа, ${\sigma }_U$ = 700 МПа, $\widetilde{{\sigma }_U}$ = 400 МПа, ${\beta }_L$ = 0.30, ${\beta }_V$ = 0.27.

Для параметров кинетического уравнения по опыту тестовых расчетов выбраны значения ${\psi }^{*}$ = 0.98, k = 0.05, $\gamma$ = 0.5, r₀ = 0.1 мм. На рис. 8, а показана геометрия образца с острой насечкой, подверженного циклическому нагружению на растяжение–сжатие. На рис. 8, b и c показаны использованные при расчете варианты измельчения сетки с характерным размером ячейки в окрестности вершины насечки 0.1 мм, рис. 8, b, и с размером 0.05 мм, рис. 8, c.

 

Рис. 8. Образец с острой насечкой без скругления – (a), сетки с характерным размером ячейки в окрестности вершины насечки 0.1 мм – (b), с размером 0.05 мм – (c).

 

Граничные условия во всех вычислительных экспериментах: циклические смещения прикладываются по продольной оси к торцам образца с амплитудой $∆$x = 26.5 мкм, касательные напряжения равны нулю. На рис. 9 показаны картины зарождения квазитрещины в вершине образца с острой насечкой, рис. 9, a, и ее развитие до длины 1 мм, рис. 9, b, и до длины 2 мм, рис. 9, c.

 

Рис. 9. Примеры расчета зарождения и развития квазитрещин. Зарождение квазитрещины – a), ее развитие до длины 1 мм – b), ее развитие до длины 2 мм – c).

 

На рис. 10, a представлены графики (пунктиром) для количества циклов N₀ до зарождения трещины на образце с остроугольной насечкой (нулевой радиус закругления) при измельчении сетки от размера r₀ до размера r₀/16, построенные по результатам расчетов с осреднением повреждаемости. На рис. 10, b эти же графики нормализованы – отнесены значению N₀ при h/r₀ = 1.

 

Рис. 10. Количество циклов N до зарождения трещины от острой насечки. Количество циклов в зависимости от размера сетки h – a), нормализованное количество циклов – b). Кривая 1 – без процедуры усреднения, кривая 2 – с усреднением.

 

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

4.3. Примеры численных расчетов зарождения и развития повреждений для различных систем дефектов в СЛП образцах. С использованием предложенной модели циклической повреждаемости и численного алгоритма были рассчитаны высокочастотные нагружения образцов с частотой нагружения 20 кГц и амплитудой 25 мкм. На рис. 11 представлена геометрия одной четверти (левой верхней) образцов с вертикальными и горизонтальными осями симметрии и системами линейных дефектов на горизонтальной оси, которые проявляются в виде концентраторов напряжений (красные точки), рис. 11, а. Дефектами в данном случае являются «непроплавленные» элементы с характерным размером 30 мкм и расстоянием между ними 100 мкм, которые имеют модуль упругости Е, составляющий одну треть от модуля основного материала образца Е₀, Е/Е₀ =0.33. На рис. 11, б представлена двухслойная система аналогичных дефектов. На этом рисунке видны две линии дефектов, проявляющих себя концентраторами напряжений при заданных параметрах нагружения (красные точки).

 

Рис. 11. Системы дефектов в образце. Концентрация напряжений на линейной системе дефектов – a), концентрация напряжений на двухслойной системе дефектов – b).

 

На рис. 12 красным цветом изображены зоны полного разрушения (квазитрещины) в момент зарождения при количестве циклов N = 1.31Е+09, рис. 12, а, и в момент, когда квазитрещина продвинулась почти на половину расстояния до оси образца при N = 1.33Е+09, рис. 12, b. Эти зоны по расположению и количеству циклов практически идентичны и для линейной, и двухслойной системы дефектов. Число циклов N =1.33Е+09 можно считать соответствующим макроразрушению образца, так как дальнейшее продвижение квазитрещины происходит практически без изменения этого значения с шагом ~10⁵–10⁶ циклов. Большее число слоев с дефектами не рассматривалось, так как более удаленные от горизонтальной оси образца дефекты практически не влияют на процесс развития повреждений, поскольку расположены в более широкой части образца с меньшим уровнем амплитуд напряжений.

 

Рис. 12. Зарождение и развитие квазитрещины на двухслойной системе дефектов. Зарождение, N =1.31Е+09 – a), развитие, N =1.33Е+09 – b).

 

Таким образом, описанная вычислительная процедура позволяет по параметрам микроструктуры СЛП материала – заданной геометрии расположения дефектов, их размеров и коэффициентов Е/Е₀, отражающих степень их непроплава-переплава, определять долговечность образцов, напечатанных методом селективного лазерного плавления. Поскольку, как показано в начальных разделах, параметры микроструктуры могут быть рассчитаны в соответствии с характеристиками мощности и скорости подвижного лазерного луча, то предложенный комплексный алгоритм позволяет оценить усталостную прочность и долговечность материала с учетом процесса селективного лазерного плавления.

4.4. Оценка количества циклов до разрушения СЛП образцов для различных систем дефектов. Систематические расчеты высокочастотного циклического нагружения образцов при заданной частоте и амплитуде позволили получить значения долговечности до разрушения в зависимости от отношения Е/Е₀ для системы дефектов, которые представлены в табл. 2.

 

Таблица 2. Результаты расчета долговечности образца для различных типов дефектов

Дефект Е/Е0	Долговечность
0.01	2.11е8
0.10	2.85е8
0.25	4.72е8
0.50	1.08е9
1.00	4.36е9
1.05	4.31е9
1.10	2.78е9
1.15	1.87е9
1.20	1.33е9

 

Значение Е/Е₀ = 1 соответствует однородному материалу без дефектов. Значения Е/Е₀ < 1 соответствуют дефектам в виде непроплавленных зон (примерно, как на рис. 7, а) с различной степенью «недоплава», а значения Е/Е₀ > 1 соответствуют дефектам в виде переплавленных зон (как на рис. 7, b). В этом случае реальные значения Е/Е₀ не превышают 1.20–1.25. Поэтому можно сказать, как это из физических соображений и ожидалось, что долговечность до разрушения материала с системами непроплавленных дефектов может уменьшаться на порядок и более (сравним значения 3.75Е+09 при Е/Е₀ =1 и 3.59Е+08 при Е/Е₀ = 0.01). Это падение близко к наблюдаемому при усталостных высокочастотных испытаниях алюминиевых сплавов [6, 7]. При появлении системы переплавленных дефектов долговечность материала до разрушения уменьшается в худшем случае в 1.5–2 раза (3.75Е+09 при Е/Е₀ =1, 2.62Е+09 при Е/Е₀ = 1.20 и 1.93Е+09 при Е/Е₀ = 1.25).

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

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

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

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

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

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

Исследование выполнено в рамках проекта РНФ № 23-19-00640.

About the authors

I. S. Nikitin

Institute for Computer Aided Design of the RAS

Author for correspondence.
Email: i_nikitin@list.ru
Russian Federation, Moscow

N. G. Burago

Institute for Computer Aided Design of the RAS; Ishlinsky Institute for Problems in Mechanics RAS

Email: i_nikitin@list.ru
Russian Federation, Moscow; Moscow

A. D. Nikitin

Institute for Computer Aided Design of the RAS

Email: i_nikitin@list.ru
Russian Federation, Moscow

B. A. Stratula

Institute for Computer Aided Design of the RAS

Email: i_nikitin@list.ru
Russian Federation, Moscow

References

Supplementary files