RESIDUAL STRESSES IN A THERMOELASTIC CYLINDER RESULTING FROM LAYER-BY-LAYER SURFACING

Full Text

Введение

Селективная лазерная плавка и спекание (SLM и SLS) - перспективные технологии изготовления трехмерных металлических деталей сложной формы с заданной неоднородностью физико-механических свойств [1–5]. Сущность технологий заключается в том, что деталь создается последовательно, шаг за шагом, за счет плавления и спекания металлических порошков в области с наперед заданной геометрической формой [6–9].

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

В современной механике сплошных сред существует ряд различных подходов к изучению феномена роста. На сегодняшний день опубликовано большое количество работ, посвященных механике растущего твердого тела. Ссылки можно найти в обзоре [10]. Работы [11; 12] посвящены развитию геометрических методов, применяемых в механике несовместных деформаций, возникающих в результате процесса наращивания. В работах [13–16] рост исследуется как непрерывный процесс осаждения деформируемых поверхностей материала на деформируемое трехмерное тело. При дополнительных предположениях о непрерывности функций, определяющих напряженно-деформированное состояние присоединяемых поверхностей материала, непрерывный процесс роста можно рассматривать как предел последовательности дискретных процессов [17;18].

Вопросам моделирования остаточных напряжений в растущих телах посвящена обширная литература [19–24]. В приближении малых деформаций, как правило, используется инкрементальный подход, который позволяет сформулировать задачу в терминах приращений полей напряжений и деформаций. При этом уравнения баланса и законы состояния оказываются подобными соответствующим уравнениям для совместных деформаций, что дает возможность применять хорошо разработанные методы решения краевых задач механики деформируемого твердого тела [25]. Для конечных деформаций этот подход не дает никакого выигрыша перед непосредственным описанием несовместных деформаций, и в этом случае используется геометрический формализм: обратная деформация, возвращающая напряженное тело в состояние, свободное от напряжений, определяется как вложение в пространство с более общей, неевклидовой геометрией [26]. В настоящей статье рассматриваются малые термоупругие деформации, возникающие в процессе послойного наращивания, и для их описания достаточно сформулировать и решить последовательность связанных динамических задач термоупругости для каждого элементарного слоя, используя в качестве начальных данных значения полей в конце временного интервала предыдущего слоя.

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

 

1. Формализация растущего тела

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

   Вестник Самарского университета. Естественнонаучная серия. 2020. Том 26, № 3. С. 63–90

   Vestnik of Samara University. Natural Science Series, 2020, vol. 26, no. 3, pp. 63–90 65

    

   

   Рис. 1.1. Схема процесса термоупругого роста

   Fig. 1.1. Diagram of the thermoelastic growth process

    

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

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

   B0 ⊂ B1 ⊂ ... ⊂ Bk ⊂ ... ⊂ BN , (1.1)

   где B0 является начальным телом и Bk - k−я сборка, образуемая начальным телом и совокупностью k слоев. Последовательности тел (1.1) ставится в соответствие числовая последовательность моментов времени присоединения слоев

   0 < τ1 < ... < τk < ... < τN . (1.2)

   Последовательности (1.1) и (1.2) в совокупности определяют сценарий роста тела. В случае дискретного роста эволюция полей деформаций и температуры может быть определена из решений краевых задач для каждого из тел Bk . При этом начальные данные для решения на k-м шаге определяются

   финальными значениями соответствующих полей на шаге k−1 и начальными (на момент присоединения)

   данными для присоединяемого слоя. В результате получаем рекуррентную последовательность краевых

   задач

   0

    

   ∀x ∈ B0 L0u0 + B0 = 0, ∀x ∈ ∂B0 D0u0 = 0, u0|t=0 = u0,

   ...

   n

    

   n

    

   ∀x ∈ Bn Lnun + B0 = 0, ∀x ∈ ∂Bn Dnun = 0, un|t=τ

    

   = 0,

    

   0 0 ( Ln−1un−1|t=τn при x ∈ Bn−1,

   Bn = Bn−1 +

   0 при x ∈ Bn\Bn−1,

   0

    

   здесь L0, ..., Ln - дифференциальные операторы, определяемые одной и той же дифференциальной операцией (уравнениями поля), но в разных областях, D0, ..., Dn - операторы граничных условий, B0 - внешняя объемная сила и плотность источника тепла и u0, ..., un - приращения полей смещения и температуры относительно начала шага, u0 - начальные данные для первого шага. Последовательность

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

   Лычев С.А., Фекри Монтасер. Остаточные напряжения в термоупругом цилиндре...

   66Lychev S.A., Fekry M. Residual stresses in a thermoelastic cylinder resulting from layer-by-layer surfacing

    

2. Начально-краевая задача для одного шага

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

   Всюду далее будем полагать, что область, занимаемая телом, представляет собой круговой цилиндр

   Ck = Dk × (0, Lˇ),

   где Dk - открытый круг, радиус

   Rˇk которого зависит от номера шага, Lˇ

   - фиксированное

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

   Rˇ0 < Rˇ1 < ... < Rˇk < ... < RˇN .

   В этой связи для построения решений удобно использовать цилиндрические координаты (rˇ, θˇ, zˇ), которые связаны с декартовыми координатами (xˇ, yˇ, zˇ) следующим образом:

   

   rˇ = ,jxˇ2 + yˇ2,

   

   θˇ = arctan yˇ ,

   xˇ

    

   zˇ = zˇ.

   При этом элементы локальных ортонормированных (физических) базисов (er, eθ, ez ) связаны c декартовым базисом (i, j, k) соотношениями

   er = cos θˇ i + sin θˇ j, eθ = − sin θˇ i + cos θˇ j, ez = k.

    

   1. Задача теплопроводности

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

      Λ ˇ 2 ˇ ˇ

      ∇ Θk − ρκ∂tˇΘk + ̟ˇ = 0,

      Θˇ k |tˇk =0 = Tˇ ,

       

      0

       

      ∂zˇΘˇ k |zˇ=0,zˇ=Lˇ = 0, ∂rˇΘˇ k |rˇ=Rˇk = 0, k

      (2.1)

      ∇

       

      где ˇ 2 - осесимметричный оператор Лапласа относительно размерных переменных, т. е.

      ˇ 2 ∂2

      1 ∂ ∂2

      

      

      ∇

       

      = + + .

      ∂rˇ2 rˇ ∂rˇ ∂zˇ2

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

      В уравнениях (2.1)

      Θˇ k =

      k

       

      Tˇk − Tˇ0

      - избыточная температура относительно референциальной

      абсолютной температуры

      k

       

      Tˇ0, ρ - плотность массы, Λ - коэффициент теплопроводности, κ -

      теплоемкость на единицу массы при постоянной деформации, ̟ˇ Введем безразмерные переменные следующим образом:

      - удельная мощность источников тепла.

      rˇ

      r =

      Rˇ0

      zˇ

      , z =

      Rˇ0

      Θˇ k

      

      k

       

      , Θ =

      Tˇm

      Λ

      , t = tˇ, (2.2)

      0

       

      κρRˇ2

      где

      Rˇ0 - радиус начального цилиндра, к которому будут присоединяться слои,

      Tˇm - температура

      плавления.

      В безразмерных переменных (2.2) краевая задача (2.1) принимает вид:

       

      ∇

       

      Здесь 2

      2

       

      ∇ Θk − ∂tΘk + ̟ = 0,

      ∂z Θk |z=0,z=L = 0, ∂r Θk |r=R k t =0 0

      k = 0, Θ | k = Tk .

      - осесимметричный оператор Лапласа относительно безразмерных переменных

       

      (2.3)

      2 ∂2

      1 ∂ ∂2

      

      

      ∇

       

      = + + ,

      ∂r2 r ∂r ∂z2

      и введены следующие безразмерные величины:

      Rˇ2

       

      Rˇk Lˇ

      

      ̟= 0

      ΛTm

      ̟ˇ , Rk =

      

      Rˇ0

      , L =

      

      .

      Rˇ0

      Вестник Самарского университета. Естественнонаучная серия. 2020. Том 26, № 3. С. 63–90

      Vestnik of Samara University. Natural Science Series, 2020, vol. 26, no. 3, pp. 63–90 67

       

      Решение начально-краевой задачи (2.3) будем искать в форме разложения

      ∞

      Θk (r, z, t) = ) θk,n(r, z)ϕk,n(t), (2.4)

      n=1

      }n=1

       

      где {θk,n ∞ - система собственных функций дифференциального оператора

      ∇

       

      Lk[u] := 2u

      с областью определения

      2

       

      ALk = {u|u ∈ L (Ck ) ∩ C

      2

      2

       

      (Ck ) ∧ n · ∇u|Γk,1 = 0 ∧ u|Γk,2 = 0 ∧ u|r=0 = O(1)� .

      В этом соотношении L (Ck ) - гильбертово пространство функций, определенных в цилиндрической

      области Ck и интегрируемых с квадратом относительно нормы I · Ik, порождаемой скалярным произведением (·, ·)k вида

      L 2π Rk

      k

       

      IuI2 := (u, u)k =

      r r r

       

      Ck

      u2 dv = r r

      0 0

      r u2 r dr dϕ dz, (2.5)

      0

      C2(Ck ) - множество всех дважды непрерывно дифференцируемых функций в области Ck , O(1) - порядковый символ Э. Ландау, означающий ограниченность функций на оси этой области, Γk,1 - боковая часть границы, а Γk,2 - основания цилиндрической области Ck , Γk,1 ∪ Γk,1 = ∂Ck, n - внешняя единичная нормаль к границе Γk,1.

      Собственные функции находятся как решения задачи Штурма–Лиувилля:

      2 ∂

      

      ∇ θk,n = λk,nθk,n,

      θk,n

      ∂r r=Rk

      = 0, θk,n

      z=0,L

      = 0.

      Нетривиальные решения этой краевой задачи могут быть представлены в виде

      

      πL

       

      R √

       

      θk,0,0 = 1 ,

      k

      

      θk,0,q =

      q

       

      J0 (r/Rkj1 )

      , q = 1, 2, 3, . . . ,

      q )

       

      J0 (j1

      

      √

       

      θk,p,0 = cos(pπz/L) , p = 1, 2, 3, . . . ,

      Rk πL/2

      (2.6)

      √

       

      θk,p,q = cos(pπz/L)

      Rk πL/2

      

      q

       

      1

       

      J0 (r/Rkj1 ) J0 (jq )

      , p = 1, 2, 3, . . . , q = 1, 2, 3, . . . .

      Им соответствуют собственные значения

      λk,0,0 = 0,

      λk,p,0 = (πp/L)2, p = 1, 2, 3, . . . ,

       

      (2.7)

      λk,p,q = (j1 2 2 2

      q ) /Rk + (πp/L) , p = 1, 2, 3, . . . , q = 1, 2, 3, . . . .

      В формулах (2.6) и (2.7) J0(x) - функция Бесселя первого рода нулевого порядка, которые для достаточно больших значений x могут быть вычислены по простым асимптотическим формулам:

      ∞

      J0(x) = )

      (−1)m x 2m

      

      ≈

      (8x − 1) cos x + (8x + 1) sin x

      √ ,

      q �∞

      m=0

      (m!)2 2

      8 πx3/2

      а {j1

      q=1

      - последовательность нулей функции Бесселя первого рода первого порядка J1(x), которые

      определяются следующими соотношениями:

      ∞

      J1(x) = )

      

      (−1)m

      x 2m+1

      ≈

      (6 − 16x) cos x + (6 + 16x) sin x

      √ .

      m=0

      m!(m + 1)! 2

      16 πx3/2

      Приведенные выше асимптотические представления для функций Бесселя имеют погрешность менее

      0.1 % при x > 30. Использование этих представлений вместо непосредственного вычисления рядов для больших значений x позволяет построить эффективные алгоритмы вычисления частичных сумм (2.4), представляющих распределения температуры и связанных с ней полей.

      Значения нулей функции J1(x) могут быть найдены численно из решения трансцендентного уравнения

      J1(x) = 0.

      Лычев С.А., Фекри Монтасер. Остаточные напряжения в термоупругом цилиндре...

      68Lychev S.A., Fekry M. Residual stresses in a thermoelastic cylinder resulting from layer-by-layer surfacing

       

      q

       

      Процедура вычисления последовательности значений j1

      без пропусков и с любой желаемой точностью

      может быть эффективно реализована, например, методом Ньютона, если в качестве начальных приближений использовать асимптотические выражения

      −1

      j1 1 J0(x) .

      q ≈ x +

      x − J1(x)

       

       

      x=qπ+π/4

      L

       

      Поскольку операторы Lk являются самосопряженными относительно метрики (2.5), то совокупность их собственных функций ортогональна и образует полную систему. Легко проверить, что нормы всех элементов системы (2.6) равны единице. Таким образом, система функций (2.6) образует ортонормиованный базис гильбертова пространства 2(Ck ) и может быть использована для построения решения начально-краевой задачи, которое будет сходиться в среднеквадратичном к заданным начальным данным и к правой части для любых фиксированных значений времени t. Несколько собственных функций и их производных по r и z изображены на рис. 2.1. Легко видеть, что их значения обращаются в ноль на основаниях области Ck , а на ее боковой поверхности обращается в ноль производная по нормали.

       

      	θ0,1,1	θ0,2,3	θ0,3,5
      θk,p,q
      ∂rθk,p,q
      ∂zθk,p,q

       

      Рис. 2.1. Собственные функции и их производные

      Fig. 2.1. Eigen functions and their derivatives

       

      Располагая явными выражениями для ортонормированной системы собственных функций, запишем разложение (2.4) в более подробном виде:

      1 ∞ J0 r 1

      Θk(r, z, t) =

      

      

      √ ϕ0,0(t) + )

      Rk jq

      

      √

      ϕ0,q (t) +

      Rk πL

      q=1 Rk

      q

       

      πLJ0(j1)

      

      I A

       

      ∞ cos pπz

       

      ∞ ∞ cos pπz

      J0 r 1

      

      + ) L

      

      ϕp,0(t) + ) ) L

      

      Rk jq

      ϕp,q (t) .

      p=1 Rk ,jπL/2

      q=1 p=1 Rk ,jπL/2

      1

       

      J0(jq )

      (2.8)

      

      II A

      

      

      I I I A

      Вестник Самарского университета. Естественнонаучная серия. 2020. Том 26, № 3. С. 63–90

      Vestnik of Samara University. Natural Science Series, 2020, vol. 26, no. 3, pp. 63–90 69

       

      Группа слагаемых, обозначенная символом ” I”, определяет осесимметричное распределение температуры, которое далее мы будем называть ’одномерной задачей’, группа слагаемых ” II” в отдельности от остальных также характеризует одномерное распределение, но, в отличие от ” I”, постоянное вдоль координаты r и переменное вдоль z, а группа слагаемых ” III” описывает двумерное распределение общего вида.

      Двухиндексная нумерация собственных значений и собственных функций, возникшая в результате разделения переменных r и z, неудобна в дальнейших вычислениях, и её предлагается заменить нумерацией с одним индексом, посредством которой собственные функции и собственные значения оказываются линейно упорядоченными

       

      n → (p, q) : (p, q) =

      /

      

      2

       

      n − 1 − 1 I

      

      

      √

       

      8n 7+1 I ,

       

       

       

      8n−7−1 I I√ −

      2 2

       

      √ √ \

      n − 1 I 8n−7+1 I I 8n−7+3 I

      , (2.9)

      2 2 2

       

      где ⌊m⌋ - нижняя функция m, которая возвращает наибольшее целое число, меньшее или равное m. Это отображение определяет пары (p, q) в порядке ”диагонального заполнения”. Его иллюстрирует следующая таблица, в ячейках которой указаны значения n, а номерам столбцов и строк отвечают

      значения (p, q).

       

      ❍ ❍❍ p q ❍❍	0	1	2	3	. . .
      0	1	3	6	10	. . .
      1	2	5	9	14	. . .
      2	4	8	13	19	. . .
      3	7	12	18	25	. . .
      ...	...	...	...	...	. . .

      Пересчет индексов (p, q) через общий индекс n позволяет представить последовательности собственных функций и собственных значений в виде линейно упорядоченных множеств θk,n, λk,n, n = 0, 1, 2, . . .:

      θk,1 = θk,0,0 λk,1 = λk,0,0

      θk,2 = θk,1,0 λk,2 = λk,1,0

      θk,3 = θk,0,1 λk,3 = λk,0,1

      . . . . . . .

      При построении решения в замкнутом виде желательно все элементы решения выразить в аналитической форме. В рамках настоящей работы будем полагать, что начальные данные и правые части уравнений (для фиксированных моментов времени t) представлены полиномами от переменных r, z в квадрате (a, b) × (c, d) ⊂ (0, Rk) × (0, L). Тогда коэффициенты Фурье этих функций относительно базисных элементов θk,n могут быть представлены линейными комбинациями интегралов

      2π L Rk d b

      r r r

      Ik,p,q (α, β; a, b; c, d) :=

      0 0 0

      r

      rαzβ Ξr (r; a, b)Ξz (z; c, d)θk,p,q (r, z) r dr dz dϕ = 2π

      c

      r

      rα+1zβθk,p,q (r, z) dr dz,

      a

      где Ξr (r; a, b) - характеристическая функция интервала (a, b) ⊂ (0, Rk ), Ξz (z; c, d) - характеристическая функция интервала (c, d) ⊂ (0, L). Эти интегралы могут быть вычислены в замкнутой форме:

      

      2√π(bα+2 −aα+2 )(dβ+1−cβ+1 )

      

      Ik,0,0(α, β; a, b; c, d) =

      √ ,

      L(α+2)(β+1)Rk

      √π Γ( α +1)(dβ+1−cβ+1 ) α α

       

      b j1 2

      

      

      

      Ik,0,q (α, β; a, b; c, d) = √ 2

      × bα+2 1F2

      + 1; 1,

      

      + 2; − 1 q 1 −

      q

       

      L(β+1)Rk J0 (j1 )

       

       

      a j1

      2 2

      2

      2Rk

      −ak+2 1F2 α

      α 1 q 1

      2 + 1; 1, 2 + 2; −

      2Rk ,

      

      √2π(bα+2 −aα+2 ) β+1

      ( icpπ

       

      icpπ

      

      L

       

      Ik,p,0(α, β; a, b; c, d) =

      (α+2)√ Rk c

      L

       

      E−β

      −

       

      − ) + E−β ( L )l

       

      −dβ+1

       

      E−β

      idpπ

      

      — L

       

      + E−β

      idpπ

      

      ,

      L

      Лычев С.А., Фекри Монтасер. Остаточные напряжения в термоупругом цилиндре...

      70Lychev S.A., Fekry M. Residual stresses in a thermoelastic cylinder resulting from layer-by-layer surfacing

       

      √

       

       

      Ik,p,q (α, β; a, b; c, d) = √ 2π

      cβ+1 (E

      (− icpπ ) + E

      ( icpπ )) −

      

      q

       

      2 LRk J0 (j1 )

      

      −β L

      

      −β L

      −dβ+1 E

      idpπ + E

       

      idpπ

      ×

      −β − L

       

      −β L

      bj1 2

       

      (2.10)

      

      ×Γ ( α + 1)

      bα+2 1F2

      

      

      

      α + 1; 1, α + 2; − 1 q 1 −

      2 2 2

      2Rk

       

      a j1

       

      2

      −aα+2 1F2 α

      α 1 q 1

      2 + 1; 1, 2 + 2; −

      2Rk .

       

      Здесь 1F2(. . .) - гипергеометрическая функция, Eα(z) - интегральная показательная функция ((α)β = α(α + 1) . . . (α + β − 1)):

       

      ∞

      1F2(x; y1, y2; z) = )

      

      (x)k zk

       

      , Eα(z) =

      

      r

       

      ∞ e−zt

       

      dt.

      k=0

      (y1)k (y2)k k! tn

      1

      Таким образом, функции, представляющие собой многочлены от букв r, z в области, представляющей собой полый цилиндр, содержащийся в Ck и продолженные нулевыми значениями за пределами этого полого цилиндра, могут быть представлены в форме разложений, все слагаемые которых определяются в замкнутой форме. Вместе с тем в практических вычислениях реализуются не полные вычисления, а соответствующие частичные суммы. Обсуждение вопроса об аналитическом выражении остаточного члена для таких частичных сумм выходит за пределы настоящей работы, однако для его оценки можно воспользоваться вычислительным анализом "практической сходимости" на тестовых примерах

      разложений. В качестве такого примера приведем "неудобную" функцию, имеющую конечные разрывы по всей границе области её ненулевых значений:

      ( z, (r, z) ∈ (a, b) × (c, d)

      θtest(r, ϕ, z) =

      0, иначе.

      Коэффициенты разложения для этой функции определяются зависимостями (2.10). Для Rk = 1, L = 3, a = 1/3, b = 2/3, c = 0, d = 2 распределение коэффициентов Фурье представлено на рис. 2.2. Линейный порядок (нумерация) членов разложения определяется зависимостями (2.9). Поверхности, построенные по частичным суммам различных порядков F(n) (число слагаемых указано на рисунке), показаны на рис. 2.3 в левой колонке. В правой колонке показаны разности значений частичных сумм и оригинальной функции. Видно, что с увеличением порядка частичной суммы эти разности уменьшаются в регулярных точках, а в точках разрыва имеются неустранимые погрешности, аналогичные известному в теории рядов Фурье “эффекту Гиббса“.

       

      

      а б

      Рис. 2.2. Распределение коэффициентов Фурье

      Fig. 2.2. Distribution of Fourier coefficients

       

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

      Вестник Самарского университета. Естественнонаучная серия. 2020. Том 26, № 3. С. 63–90

      Vestnik of Samara University. Natural Science Series, 2020, vol. 26, no. 3, pp. 63–90 71

       

      	F(n)	θtest − F(n)
      n = 54 собственных функций
      n = 216 собственных функций
      n = 900 собственных функций

       

      Рис. 2.3. Частичные суммы разного числа собственных функций

      Fig. 2.3. Partial sums of different numbers of eigen functions

       

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

      При известных собственных функциях θk,n(r, z) координатные функции ϕk,n(t) могут быть найдены из решения несвязанной последовательности задач Коши

       

       

      λϕk,n − ϕ˙ k,n = Ωk,n, ϕk,n

       

       

      t=t0

      = Φk,n,

       

      где Ωk,n = Ωk,n(t) определяется проекциями заданной мощности источников тепла ωk (r, z, t) на собственные функции θk,n(r, z), т. е.

       

      L Rk

      r r

      Ωk,n(t) = 2π

      0 0

      ωk(r, z, t) θk,n(r, z) r dr dz,

      Лычев С.А., Фекри Монтасер. Остаточные напряжения в термоупругом цилиндре...

      72Lychev S.A., Fekry M. Residual stresses in a thermoelastic cylinder resulting from layer-by-layer surfacing

       

      а Φk,n - проекциями начальных данных:

       

      L Rk

      r r

      Φk,n = 2π

      0 0

      k

       

      Θ0 (r, z, t) θk,n(r, z) r dr dz.

      Решение этих задач может быть представлено в виде

      t

      r

      ϕk,n(t) = Φk,n eλk,n (t−t0 ) +

       

      Ωk,n(τ ) eλk,n (t−τ ) dτ. (2.11)

      t0

      В случае постоянного распределения Ωk,n, т. е. при Ωk,n(τ ) = Ωk,n = const,

      ϕk,n(t) = Φk,neλk,n(t−tn ) + Ωk,n(1 − et−tn ).

      Таким образом, на отдельно взятом шаге распределение температуры в теле полностью определяется разложением (2.8). Для иллюстрации эффективности такого представления рассмотрим пример разложения решения однородного уравнения (2.5) (при ω = 0) и для начального распределения температуры, которое определяется “bump“-функцией, т. е. гладкой функцией переменных r, z, обращающейся в ноль на границе области, занимаемой цилиндром (r, z) ∈ (0, 1) × (0, 3):

      Θ0(r, z) = 4096r3(1 − 3r + 3r2 − r3)(z/3)3(1 − z + z2/3 − z3/27).

      Поверхности, представляющие решения в плоскости параметров (r, z) для следующих значений безразмерного времени t = 0, 0.015, и 0.03, показаны на рис. 2.4, а, б, в соответственно. Ниже, на рис. 2.4, г, д, е для наглядности показаны распределения полей температуры в трехмерном пространстве.

       

      

      

      

      а б в

      

      г д е

      Рис. 2.4. Поля температуры в различные моменты времени. Решение тестовой задачи.

      Температура представлена в безразмерной форме и изменяется от 0 до 1. Шаг изолиний - 0.1 Fig. 2.4. Temperature fields at different times. Test problem solution. The temperature is presented in dimensionless form and varies from 0 to 1. The isoline step is 0.1

       

      Для вычисления распределения температуры на всех шагах процесса k = 1, 2, 3, . . . , N требуются еще формулы для вычисления коэффициентов Фурье начальных данных, взятых с финального момента

      Вестник Самарского университета. Естественнонаучная серия. 2020. Том 26, № 3. С. 63–90

      Vestnik of Samara University. Natural Science Series, 2020, vol. 26, no. 3, pp. 63–90 73

       

      Θ˜ 0

       

      предыдущего шага, Θk−1(r, z, tk), и продолженных распределением температуры в присоединенном слое,

      k(r, z), т. е.

      Ck,n =

      r r r

      Θk−1(r, z, tk)θk,n(r, z) dv +

      r r r

      Θ˜ 0

       

      k(r, z) θk,n(r, z) dv =

       

      L Rk−1

      Ck−1

      Ck \Ck−1

       

      L Rk

      r r

      = 2π

      ∞

      r

       

      r

       

      ) θk−1,s(r, z) ϕk−1,s(tk ) θk,n(r, z) r dr dz + 2π

      Θ˜ 0

       

      k(r, z) θk,n(r, z) r dr dz.

      0 0 s=1

      0 Rk−1

      Некоторая техническая сложность состоит в том, что распределение температуры в финальный момент (k−1)-го шага представлено в форме разложений по собственным функциям оператора Lk−1, в то время как для решения уравнения (2.11) требуется разложение по собственным функциям оператора Lk. Оба набора собственных функций выражены через тригонометрические функции и функции Бесселя, однако последовательности собственных значений различаются. Эта сложность может быть преодолена, если воспользоваться формулами:

      Rk−1

      3 1

      r r

      1 J

      j

       

      r 1 r dr = Rk−1Rkjq .

      Jk,p,q =

      

      J0 jp

      Rk−1

      

      

      0 Rk q

      R2 (j1 2

      1 (j1)2

      0 k−1

      q ) − Rk p

      Таким образом, если распределение температуры в каждом присоединяемом слое в момент присоединения задано полиномиальным распределением

       

      Θ˜ 0

       

      k(r, z) =

      ),

      α,β∈N

      A˜k,α,β rαzβ, (r, z) ∈ (a′, b′) × (c′, d′)

       

      (2.12)

       

      где

      0, иначе,

      A˜k,α,β - заданные числа, а интервалы (a′, b′) ⊂ (Rk−1 , Rk), (c′, d′) ⊂ (0, L) определяют область

      внутри присоединяемого слоя, то

      ∞ L Jk,q(n),q(s) ϕk−1 (tk ) ( 0, p(n) I= p(s)

      Ck,n = 2π )

      1/2, p(n) = p(s) I= 0

      +2π )

      Ak,α,β Ik,p(n),q(n)(α, β; a′, b′, c′, d′).

      s=1 Rk Rk−1J0

      j

       

      1

      q(n)

      j

       

      J0

       

      1

      q(s)

      1, p(n) = p(s) = 0

      α,β∈N

       

      (2.13)

      Окончательно приходим к следующей рекуррентной последовательности решений, которые определяют распределение температуры в дискретно растущем цилиндре в течение всего процесса роста. Полагаем, что температура начального тела - цилиндра радиуса R0 - распределена по следующему закону:

       

      Θ0

       

      0(r, z) =

      ),

      α,β∈N

      A0,α,βrαzβ, (r, z) ∈ (a0, b0) × (c0, d0)

      0, иначе.

      Здесь (a0, b0) ⊂ (0, R0), (c0, d0) ⊂ (0, L) - интервалы, определяющие часть начального тела C0, избыточная температура которой отлична от нуля. Распределение избыточной температуры в слоях задается зависимостями (2.12). Мощность источников тепла распределена в начальном теле C0 следующим образом:

      ),

      B0,α,βrαzβψ0,α.β (t), (r, z) ∈ (a′′, b′′) × (c′′, d′′)

      ω0(r, z, t) =

      0 0 0 0

      α,β∈N

      0, иначе,

      (2.14)

      где (a′′, b′′) ⊂ (0, R0), (c′′, d′′) ⊂ (0, L) - интервалы, определяющие часть начального тела C0, мощность

      0 0 0 0

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

      телом, увеличивается, и для описания распределения мощности источников тепла на каждом шаге задается зависимость, аналогичная зависимости (2.14), однако уже в расширенной области и, возможно, с другими коэффициентами разложения (k = 1, 2, . . . , N )

      ),

      Bk,α,βrαzβψ0,α.β (t), (r, z) ∈ (a′′, b′′) × (c′′, d′′)

      ωk(r, z, t) =

      k k k k

      α,β∈N

      0, иначе,

      причем интервал (a′′, b′′) может быть частью интервала более широкого, чем в начале процесса:

      k k

      k, bk ) ⊂ (0, Rk ), а интервал (ck, dk ) по-прежнему вложен в интервал (0, L).

      (a′′ ′′

      ′′ ′′

      Решение на каждом шаге представляется разложением (2.4), в котором базисные функции θ0,n(r, t) определяются по (2.7), а координатные функции ϕ0,n(r, t) находятся рекуррентным вычислением. Базой рекурсии являются координатные функции, найденные для начального тела

      ϕ0,n(t) = )

      Aα,β I0,n(α, β; a0, b0; c0, d0) eλ0,nt + )

      t

      r

      B0,α,β I0,n(α, β; a′′, b′′; c′′, d′′)

       

      ψ0,α,β (τ )eλ0,n (t−τ ) dτ,

      α,β∈N

      α,β∈N

      0 0 0 0

      0

      Лычев С.А., Фекри Монтасер. Остаточные напряжения в термоупругом цилиндре...

      74Lychev S.A., Fekry M. Residual stresses in a thermoelastic cylinder resulting from layer-by-layer surfacing

       

      а координатные функции для последующих шагов находятся по формуле

      t

      ϕk,n(t) = Ck,n eλk,n (t−ts ) + )

      r

      Bk,α,β I0,n(α, β; a′′ , b′′; c′′ , d′′)

      ψk,γ,δ (τ )eλk,n (t−tk −τ ) dτ,

      γ,δ∈N

      n n n n

      tk

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

       

   2. Характерное время процесса теплопроводности

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

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

      Характерное время технологического процесса наращивания может быть определено как интервал между моментами присоединения смежных слоев. Характерное время процесса теплопроводности предлагается определить как время, за которое температура на границе роста системы ”цилиндр-слой” уменьшится в e раз относительно температуры в момент присоединения слоя. Для определения этого времени достаточно воспользоваться подстроенным выше решением для случая N = 1 (т. е. системы "начальный цилиндр-слой"). Тогда характерное время процесса теплопроводности tth будет определяться как

      tth = t∗ − t1,

      где t1 - момент присоединения слоя, t∗ - наименьший корень уравнения

      Θ1(R1, z0, t∗) = Θ1(R1, z0, t1)/e,

      в котором z0 определяет некоторую фиксированную координату по z, задающую плоскость измерения. В расчетах принималось z0 = L/2.

      Значение характерного времени процесса теплопроводности, помимо физико-механических свойств материала, зависит от соотношения толщины слоя и радиуса начального тела. В настоящей работе рассматривались слои, толщина которых составляет 1 % от радиуса начального цилиндра, который, в свою очередь, принимался равным 10, 1 и 0.1 мм. В безразмерных переменных характерное время не зависит от свойств материала и размера тела, и определяется только отношением толщины слоя к радиусу цилиндра. Для указанного отношения оно составляет 0.000222287. Размерные величины характерных времен приведены для титана и меди в таблице.

       

       

      Радиус, мм	10	1	0.1
      Характерное время для титана, с Характерное время для меди, с	0.310679 0.0213001	0.00310679 0.000213001	0.0000310679 2.13001 · 10−6

       

      Размерные времена процессов для различных материалов и размеров тел Dimensional times of processes for various materials and sizes of bodies

      Таблица

       

      Table

       

      Используя характерное время процесса теплопроводности как главную характеристику темпа связанного процесса теплопроводности и роста, можно разделить технологические процессы на три класса:

      • Быстрые (интервал между моментами присоединения существенно меньше характерного времени).

      • Соразмерные (шаг процесса одного порядка с характерным временем).

      • Медленные (интервал между моментами присоединения существенно больше характерного времени).

        Ниже приведены результаты вычислительного моделирования, выполненные на основе построенного выше рекурсивного решения. Безразмерный радиус принимался равным 1, высота 3, а толщины слоев - 0.01, т. е. 1 % от начального радиуса. В расчетах учитывались 15 слоев, а интервалы времени между

        Вестник Самарского университета. Естественнонаучная серия. 2020. Том 26, № 3. С. 63–90

        Vestnik of Samara University. Natural Science Series, 2020, vol. 26, no. 3, pp. 63–90 75

         

        моментами присоединения соседних слоев выбирались так, чтобы аддитивный процесс соответствовал одному из приведенных выше типов. В первой серии расчетов, которая соответствует “быстрому процессу“, этот интервал принимался равным характерному времени, уменьшенному в 100 раз, т. е. tth/100, во второй – равным характерному времени (“соразмерный процесс“), а в третьей - в 100 раз большим, т. е. 100tth (“медленный процесс“). Момент присоединения первого слоя полагался равным утроенному типовому шагу. Схематичная диаграмма процесса показана на рис. 2.5.

         

        

         

        

         

        Рис. 2.5. Схематичная диаграмма процесса роста

        Fig. 2.5. Schematic diagram of the growth process

         

        В частичных суммах разложений учитывались 300 слагаемых. Температура присоединяемых слоев в момент присоединения предполагалась постоянной во всех точках слоя и равной 1, т. е. безразмерной температуре плавления. Безразмерная температура начального цилиндра в начальный момент времени считалась равной 0, т. е. лабораторной температуре. Мощность источников тепла, подводимого для дополнительного прогрева, задавалась распределением

         0, r < Rk − δ

        

        

        ωk = 

        

        Pk 0.0336785 −

        5.1608240ρ + 298.23072ρ2

        − 6256.5729ρ3

         

        , δ < r < Rk.

         

         

        −

         

        

         

         +64398.198ρ4 310609.02ρ5 + 623226.03ρ6

        ρ=r−Rk +δ

        В этом распределении δ характеризует глубину зоны индукционного прогрева, которая зависит от частоты возбуждающего тока, Pk - мощность прогрева на шаге k. Значения при k = 0 соответствуют прогреву начального тела. Задавая различные зависимости для n 1→ Pn, можно

        реализовывать различные сценарии коррекции аддитивного процесса. В расчетах использовались

        следующие зависимости для Pn.

        Для “быстрого процесса“ полагалось

         

        Pk =

        ( 1.2 · 104, k = 0,

        0, k > 0,

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

        температуры плавления и не возникал перегрев её окрестности, т. е. температура внутри области не превышала температуру плавления (в безразмерных переменных - 1).

        Для “соразмерного процесса“ использовался “закон обратных квадратов“:

        ( 1.8 · 102/k2, k < 15,

        Pk =

        0, k = 15.

        В момент присоединения последнего слоя индукционный нагрев прекращается. Для “медленного процесса“ использовалась обратная пропорциональность:

        ( 8/k2, k < 15,

        Pk =

        0, k = 15.

        Столбчатые диаграммы, характеризующие распределения Pn для процессов различного типа, показаны на рис. 2.6. Поверхности, характеризующие изменение температурного профиля в течение аддитивного процесса для трех типов процессов, показаны на рис. 2.7. В левой колонке показаны распределения

        Лычев С.А., Фекри Монтасер. Остаточные напряжения в термоупругом цилиндре...

        

        

        

        76Lychev S.A., Fekry M. Residual stresses in a thermoelastic cylinder resulting from layer-by-layer surfacing

         

        Быстрый процесс Соразмерный процесс Медленный процесс Рис. 2.6. Распределение источника тепла

        Fig. 2.6. Distribution of the heat source

         

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

         

        	Без нагрева	С нагревом
        Быстрый процесс
        Соразмерный процесс
        Медленный процесс

         

        Рис. 2.7. Распределение температуры

        Fig. 2.7. Temperature distribution

        Вестник Самарского университета. Естественнонаучная серия. 2020. Том 26, № 3. С. 63–90

        Vestnik of Samara University. Natural Science Series, 2020, vol. 26, no. 3, pp. 63–90 77

         

        

        Быстрый процесс

         

        Без нагрева С нагревом

         

        Медленный процесс

         

        Соразмерный процесс

         

        Рис. 2.8. Распределение температуры на движущейся границе

        Fig. 2.8. Temperature distribution on the moving boundary

         

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

         

   3. Остаточные напряжения

Для анализа напряженного состояния, вызванного неравномерным нагревом и наращиванием,

рассмотрим задачу термоупругости, сформулированную относительно приращений перемещений

uˇk на

каждом шаге. В приближении малых деформаций квазистатические уравнения могут быть записаны в виде [31; 32]:

µ ˇ 2

∇ ˇ ˇ ˇ ˇ

∇ · uˇk − γ∇Ξk + Xk = 0, (2.15)

где Ξˇk = Ξˇk (tˇ) = Θˇ k(tˇ)−Θˇ k (tˇk ) – приращение температуры на k-м шаге, γ = (3λ+2µ)α, α - коэффициент линейного теплового расширения, λ и µ - упругие модули Ламе, ρ, как и выше, плотность массы, а Xk - приращение плотности объемных сил на k-м шаге. Уравнения сформулированы относительно размерного приращения перемещений на k-м шаге, которое, согласно договоренности, принятой выше, обозначено верхним символом ‘ ˇ ‘.

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

σˇ k = [λ ˇ · uˇk − γΞˇk ]I + µ[ ˇ uˇk + ( ˇ uˇk )T ].

∇ ∇ ∇

Для постановки краевой задачи следует указать краевые и начальные условия. Они могут быть записаны следующим образом:

 

n · σˇ k

 

 

Γ 1

 

= 0, n ⊗ n · uˇk + n · σˇ k · (I − n ⊗ n)

 

 

Γk,2

= 0,

uˇk|t=0 = 0.

k,

 

Вводя безразмерные независимые и зависимые переменные, по формулам (2.2) и в дополнение к ним, полагая uk = uˇk /Rˇ0, сформулируем начально-краевую задачу в безразмерной форме:

2 uk

∂ 1 ∂

∂wk ∂

∇ uk − r2 + K ∂r

(r uk ) +

r ∂r ∂r

− γ ∂r Ξk = 0, (2.16)

Лычев С.А., Фекри Монтасер. Остаточные напряжения в термоупругом цилиндре...

78Lychev S.A., Fekry M. Residual stresses in a thermoelastic cylinder resulting from layer-by-layer surfacing

 

	Без нагрева	С нагревом
Быстрый процесс
Соразмерный процесс
Медленный процесс

 

Рис. 2.9. Изменение температурного профиля после окончания процесса

Fig. 2.9. Changing the temperature profile after the end of the process

 

2 ∂ 1 ∂

∂wk ∂

 

где

∇ wk + K ∂z

(r uk ) +

r ∂r ∂r

− γ ∂z Ξk = 0, (2.17)

λ λ

K = µ + 1, γ = α

3 + 2

µ

, Ξk (t) = Θk(t) − Θk (tk−1).

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

 

Γ1 : ∂wk

 

∂uk

 

= 0, 2 ∂uk

 

∂uk

 

uk ∂wk

 

= γΞk,

∂r +

 

 

∂z

r=Rk

∂r + (K − 1) ( ∂r + r +

)

∂z

r=Rk

 

 

Γ2 : wk

 

 

= 0, ∂wk + ∂uk

= 0.

z=0,L

∂r ∂z z=0,L

Имея в виду условия скользящего контакта на основаниях цилиндра, решение краевой задачи будем

Вестник Самарского университета. Естественнонаучная серия. 2020. Том 26, № 3. С. 63–90

Vestnik of Samara University. Natural Science Series, 2020, vol. 26, no. 3, pp. 63–90 79

 

разыскивать в форме разложения Фурье по координате z:

∞ ∞

uk (r, z; t) = ) Uk,p(r; t) cos p˜z, wk (r, z; t) = ) Wk,p(r; t) sin p˜z,

 

где для краткости обозначено

p=0

 

p˜ =

 

∈

 

pπ , p N. L

p=1

Подстановка этих разложений в уравнения (2.17) и вычисление проекций на базисные функции

 

L

Pr f (z) := r

g(z)

0

 

f (z) cos p˜z dz g(z) sin p˜z

 

 

приводит к последовательности краевых задач, определяемых обыкновенными дифференциальными уравнениями

(K + 1) U ′′

+ 1 U ′

1

 

— Kp˜W ′ − U

K + p˜2r2 + 1

− Xk,p = 0,

k,p

W ′′

r k,p

′

k,p

′

r2 k,p (

2

)

 

p˜

 

и краевыми условиями

r

 

k,p + Kp˜Uk,p + Wk,p − (K + 1)p˜ Wk,p + K

Uk,p − Yk,p = 0

 

W ′ = 0, 2U ′

1

+ (K − 1) U ′ + U

 

+ p˜W

 

= Z , (2.18)

 

где

 

 

k,p − p˜Uk,p

r=Rk

 

L

k,p

k,p

r k,p

k,p

 

 

r=Rk

k,p

∂r

 

Xk,p = γ { ∂Ξk cos p˜z dz,

0

L

∂z

 

Yk,p = γ { ∂Ξk sin p˜z dz = γ

0

/ L \

Ξk (L) − Ξk (0) + { Ξk cos p˜dz ,

0

L

 

 

Zk,p = γ {

0

Ξk

r=Rk −0

cos p˜z dz.

Подстановка в эти выражения разложения для температуры дает

1

 

√ ∞ jq J1 ( r

1 '\

q

Xk,p = −√ γ L

1

 

), Rk j

[ϕk,p,q

(t) − ϕk,p,q

(tk−1

)] ,

k

 

(2−δp,0 )πR2 q=1

γn√πL/2

J0 (jq )

 

∞ J0 ( r 1

Yk,p = − Rk

[ϕk,p,0(t) − ϕk,p,0(tk )] + ),

q=1

Rk jq )

q )

 

J0 (j1

[ϕk,p,q (t) − ϕk,p,q (tk−1)] ,

∞

−1 )

Zk,p = γI L

), ϕk,p,q (t)−ϕk,p,q (tk .

π(2−δp,0 ) q=0 Rk

Для построения решения краевой задачи вначале найдем решения системы однородных уравнений

(K + 1) U ′′

+ 1 U ′

′

 

— Kp˜Wk,p,s − 1 U

+ p˜2r2 + 1) = 0,

k,p,s

W ′′

r k,p,s

′ ′

r2 k,p,s (K

2

 

p˜

r

 

k,p,s + Kp˜Uk,p,s + Wk,p,s − (K + 1)p˜ Wk,p,s + K

Uk,p,s = 0,

где индекс s = 1, 2, 3, 4 указывает на частное решение из фундаментальной системы решений. Эти решения могут быть представлены в виде

Uk,p,1

Wk,p,1

=

 

I1(p˜r) ,

I0(p˜r)

Uk,p,2

Wk,p,2

= p˜rI0(p˜r) − 2(1/K + 1)I1(p˜r) ,

p˜rI1(p˜r)

Uk,p,3

Wk,p,3

=

 

−K1(p˜r) ,

K0(p˜r)

Uk,p,4

Wk,p,4

= p˜rK0(p˜r) + 2(1/K + 1)K1(p˜r) .

−p˜rK1(p˜r)

Здесь Ip(. . .), Kp(. . .), p = 1, 2 - модифицированные функции Бесселя первого и второго родов. Имея фундаментальную систему решений, неоднородные уравнения могут быть решены методом вариации произвольных постоянных. В результате вычислений получим следующее выражение:

4

( U )k,p Wk,p = Ck,p,1 ( U )k,p,1 Wk,p,1 + Ck,p,2 ( U )k,p,2 Wk,p,2 + ) Ak,p,s ( U )k,p,s Wk,p,s, (2.19)

s=1

Лычев С.А., Фекри Монтасер. Остаточные напряжения в термоупругом цилиндре...

80Lychev S.A., Fekry M. Residual stresses in a thermoelastic cylinder resulting from layer-by-layer surfacing

 

где

r J

2

 

Ak,p,1 = 1 {

0

K+1

 

K0(p˜ρ) (Kp˜ρXk,p(ρ) + 2Yk,p(ρ)) − Kp˜ ρYk,p(ρ)K1(p˜ρ)

r

ρ dρ,

Ak,p,2 = − 1 { JI0(p˜ρ) (Kp˜ρXk,p(ρ) + 2Yk,p(ρ)) + Kp˜ ρYk,p(ρ)I1(p˜ρ)) ρ dρ,

r

 

2 0 K+1

Ak,p,3 = K { J 1

 

k,p 0

k,p 1

2 K+1 Y

0

r

(ρ)K (p˜ρ) − X

(ρ)K (p˜ρ)

ρ dρ,

Ak,p,4 = − 2( K 1) {

(K + 1)Xk,p(ρ)I1(p˜ρ) + Yk,p(ρ)I0(p˜ρ)} ρ dρ,

K+ {

0

а Ck,p,s находятся из решения системы линейных уравнений, получаемой при подстановке (2.19) в краевые условия (2.18), т. е.

 

Ck,p,1 =

 

 

B1 D12

 

 

B2 D22

 

 

D11 D12

 

 

D21 D22

 

 

 

 

 

, Ck,p,2 =

 

 

 

 

 

 

 

D11 B1

 

 

 

 

D12 B2

.

 

 

 

 

D11 D12

 

 

 

 

D21 D22

Элементы определителей определяются по следующим соотношениям:

D11 = B1[Uk,p,1, Wk,p,1] = 2p˜I1(p˜Rk ),

K

 

D12 = B1[Uk,p,2, Wk,p,2] = 2p˜ ( 1 + 1) I1(p˜Rk ] − p˜Rk I0(p˜Rk )l ,

Rk

 

D21 = B2[Uk,p,1, Wk,p,1] = −2p˜I0(p˜Rk ) + 2 I1(p˜Rk ),

D22 = B2[Uk,p,2, Wk,p,2] = 2 12+K(2+p˜ R ) 1

k − (p˜ + 2Kp˜)I0(p˜Rk )1 ;

2 2

k

 

K r I (p˜R )

4 4

B1 = −B1

), Ak,p,sUk,p,s, ), Ak,p,sWk,p,s ,

s=1

4

s=1

4

B2 = Zk,p − B2

), Ak,p,s Uk,p,s, ), Ak,p,sWk,p,s ,

s=1

s=1

где B1, B2 – операторы, определяемые краевыми условиями (2.18):

 

 

B1[U, W ] = W ′ − p˜U

 

 

 

r=Rk

, B1[U, W ] = 2U ′ + (K − 1)

 

U ′ +

1

U + p˜W

r

 

 

 

 

 

.

 

r=Rk

Завершающим шагом построения решения является вычисление интегралов (2.3.). Вычислить их полностью в терминах известных специальных функций не удается, поскольку интегралы от произведений функций Бесселя различных аргументов и порядков с весом ρ2 не табулированы в справочной литературе. Однако их можно привести к виду

γp˜√πL/2 ∞

 

1 j1 K

 

Kpp˜

Ak,p,1 = − Rk

),

1

 

J0 (j

q

 

√2(2

0,1(r) +

S1,0(r)+

 

+ 1

q=1

q ) −δp,0 )LRk S

 

2(K+1)

p˜2 +q˜2 [1 + rq˜J1 (q˜r) K0 (p˜r) − p˜rJ0 (q˜r) K1 (p˜r)]

[ϕk,p,q (t) − ϕk,p,q (tk )] −

 

/ 2 1

\

1−rp˜K1 (p˜r)

Kp˜

2,1

rp˜

− p˜2 + K+1 G1,3 2

1/2, 3/2, 0

[ϕk,p,0(t) − ϕk,p,0(tk )] ,

γp˜√πL/2 ∞

 

1 j1 K

Kpp˜

Ak,p,2 = Rk

),

1

 

J0 (j

q

 

√ 0,1(r) + 2(K+1)

1,0(r)+

+ r

q=1

q ) 2(2−δp,0 )LRk P P

 

p˜2 +q˜2 [p˜I1(p˜r)J0 (q˜r) + q˜I0(p˜r)J1 (q˜r)]

[ϕk,p,q (t) − ϕk,p,q (tk )] −

( rI1 (p˜r)

— p˜

+ F

 

Kp˜2 r3

6(K+1) 1 2

 

3/2; 2, 5/2;

rp˜ 2

2

[ϕk,p,0(t) − ϕk,p,0(tk )] ,

√ ∞ √

4√πR2

 

(K+1)(p˜ +q˜ )J0 (j ) +

 

Ak,p,3 = γK L

k

),

q=1

2πp˜Rk [−q˜rK0 (p˜r)J1 (q˜r)+p˜rK1 (p˜r)J0 (q˜r)−1]

2 2 1

q

q

 

2j1 [p˜2 (−r)K2 (p˜r)J1 (q˜r)+p˜q˜rK1 (p˜r)J2 (q˜r)+q˜]

q

 

+ (p˜3 +p˜q˜2 )√2−δp,0 J0 (j1 )

 

√

 

+ 2πRk (p˜rK1 (p˜r)−1)

[ϕk,p,q (t) − ϕk,p,q (tk )] +

 

(K+1)p˜

[ϕk,p,0(t) − ϕk,p,0(tk )] ,

Вестник Самарского университета. Естественнонаучная серия. 2020. Том 26, № 3. С. 63–90

Vestnik of Samara University. Natural Science Series, 2020, vol. 26, no. 3, pp. 63–90 81

 

∞

 

Ak,p,4 = ), − γK

√

Lr

 

q=1

8√π(K+1)R2 (p˜2 +q˜2 )√2−δp,0 J0 (j1 )×

k q

/ 2 2

q )

 

× πp˜Rk ,j4 − 2δp,0

4

 

r (p˜2 + q˜2) 0F˜1 ; 2; p˜ r

[ϕk,p,0(t) − ϕk,p,0(tk )] J0 (j1 +

\

+2φk,p,q [p˜I1(p˜r)J0(q˜r) + q˜I0(p˜r)J1(q˜r)] −

 

q

 

−4(K + 1)j1 [ϕk,p,q (t) − ϕk,p,q (tk )] [q˜I1(p˜r)J0 (q˜r) − p˜I0(p˜r)J1(q˜r)] ,

 

где q˜ =

j

 

1

q , G2,1(. . .) – функция Мейера [33], а Sα,β , Pα,β – функции переменной r, которые могут быть

Rk 1,2

заданы следующими формулами:

r r

r r

Sα,β (r) :=

0

Kα(p˜ρ)Jβ (q˜ρ)ρ2 dρ, Pα,β (r) :=

0

Iα(p˜ρ)Jβ (q˜ρ)ρ2 dρ.

В одномерном случае, т. е. в случае, когда искомое решение не зависит от z и W ≡ 0, в разложении (2.4) отлично от нуля только одно слагаемое – Uk,0, которое может быть найдено как решение дифференциального уравнения

(K + 1)

 

U

 

′′

k,0

1

+ U

 

′

r k,0

1

− r2 Uk,0 (K + 1) − Xk,0 = 0,

удовлетворяющее краевому условию

2U

 

U

 

1

′

k,0

+ (K − 1)

′

k,0

+ Uk,0

r

 

 

r=Rk

= Zk,0.

Решение сформулированной таким образом краевой задачи может быть представлено в виде

r

 

1 r

 

ρ2

Rk

r r ρ2

 

rZk,0

Uk,0 = 2(K + 1)

0

−

 

r Xk,0(ρ) dρ −

r K

0

R

 

K + 2

k

Xk,0(ρ) dρ + 2K ,

который после подстановки выражения разложения для Xk,0(ρ) и вычисления соответствующих интегралов может быть преобразован также к разложению по функциям Бесселя:

γ ∞ J1 r 1

Uk,0 = √ )

Rk jq

1 1

[ϕk,0,q (t) − ϕk,0,q (tk )] .

πL(K + 1) q=1

jq J0(jq )

Полученным таким образом приращениям смещений соответствуют приращения безразмерных напряжений (отнесенные к модулю сдвига)

2γ ∞

J1 r 1

√ )

Rk jq

[ϕk,0,q (t) − ϕk,0,q (tk )] + +

1

[ϕk,0,0(t) − ϕk,0,0(tk )] ,

σrr,k = − r(

+ 1) πL

j1 1

K q=1

q J0(jq ) Rk

γ 2 ∞

1Rk J1 r

1 r 1

 

— rj1J0 j1

σθθ,k = √ )

Rk jq

1

q Rk q

1

1

××[ϕk,0,q (t) − ϕk,0,q (tk )]+

[ϕk,0,0(t) − ϕk,0,0(tk )] ,

πL (K + 1)rRk

q=1

jq J0 (jq ) Rk

2γ ∞

J0 r 1

√ )

Rk jq

[ϕk,0,q (t)

ϕk,0,q k

1

R

 

k,0,0

k,0,0 k

σzz,k = − R (

+ 1) πL

J (j1

− (t )] + + [ϕ

k

(t) − ϕ

(t )] .

k K q=1

0 q )

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

4

 

k,0 − Uk,0)2 + r Uk,0Uk,0l .

σMises,k = − 3r2 (r U ′ ′

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

Лычев С.А., Фекри Монтасер. Остаточные напряжения в термоупругом цилиндре...

82Lychev S.A., Fekry M. Residual stresses in a thermoelastic cylinder resulting from layer-by-layer surfacing

 

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

Расчеты проводились для титана [34]:

λ = 113.8 ГПа, µ = 44 ГПа, ρ = 4510 кг/м3, α = 8.6 · 10−6 К−1,

Λ = 17 Вт/(м·К), κ = 521 Дж/(кг·К).

Распределения приращений безразмерных перемещений, изменяющихся в ходе процесса, показаны

на рис. 2.10.

 

	Без нагрева	С нагревом
Быстрый процесс
Соразмерный процесс
Медленный процесс

 

Рис. 2.10. Приращения безразмерных перемещений в ходе процесса

Fig. 2.10. Increments of dimensionless displacements during the process

 

Суммарные безразмерные смещения от места присоединения могут быть определены как

k−1

 

 

u = u˜k + ) u˜s ,

t=ts

s=1

где u˜s - приращения перемещений, продолженные нулевыми значениями:

( us, r < Rs

u˜s =

0, r > Rs,

Вестник Самарского университета. Естественнонаучная серия. 2020. Том 26, № 3. С. 63–90

Vestnik of Samara University. Natural Science Series, 2020, vol. 26, no. 3, pp. 63–90 83

 

а суммирование ведется до индекса k, который для рассматриваемого момента времени t удовлетворяет условию

 

tk <� t <

( tk+1, k < N

∞, k = N

 

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

 

	Без нагрева	С нагревом
Быстрый процесс
Соразмерный процесс
Медленный процесс

 

Рис. 2.11. Изменение суммарных безразмерных смещений (от места присоединения) в ходе процесса

Fig. 2.11. Change in the total dimensionless displacements (from the point of attachment) during the process

Лычев С.А., Фекри Монтасер. Остаточные напряжения в термоупругом цилиндре...

84Lychev S.A., Fekry M. Residual stresses in a thermoelastic cylinder resulting from layer-by-layer surfacing

 

	Без нагрева	С нагревом
Быстрый процесс
Соразмерный процесс
Медленный процесс

 

Рис. 2.12. Эволюция безразмерных радиальных компонент перемещений

Fig. 2.12. Evolution of dimensionless radial components of displacements

 

Заключение

Анализ результатов, полученных в ходе вычислений, позволяет сделать следующие выводы.

 

1. Индукционный нагрев может существенно изменить распределение остаточных напряжений и искажений.

    

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

    

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

   Вестник Самарского университета. Естественнонаучная серия. 2020. Том 26, № 3. С. 63–90

   Vestnik of Samara University. Natural Science Series, 2020, vol. 26, no. 3, pp. 63–90 85

    

   	Без нагрева	С нагревом
   Быстрый процесс
   Соразмерный процесс
   Медленный процесс

    

   Рис. 2.13. Эволюция безразмерных радиальных компонент перемещений

   

   

   

   Fig. 2.13. Evolution of dimensionless radial components of displacements

    

   

    

   

    

   

    

   

    

   

    

   Быстрый процесс Соразмерный процесс Медленный процесс Рис. 2.14. Остаточные безразмерные перемещения

   Fig. 2.14. Residual dimensionless displacements

   Лычев С.А., Фекри Монтасер. Остаточные напряжения в термоупругом цилиндре...

   

   

   

   86Lychev S.A., Fekry M. Residual stresses in a thermoelastic cylinder resulting from layer-by-layer surfacing

    

   

    

   Быстрый процесс Соразмерный процесс Медленный процесс Рис. 2.15. Остаточные безразмерные интенсивности напряжений

   Fig. 2.15. Residual dimensionless stress intensities

    

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

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

 

About the authors

S. A. Lychev

IPMech RAS

Author for correspondence.
Email: lychevsa@mail.ru
ORCID iD: 0000-0001-7590-1389

Doctor of Physical and Mathematical Sciences, associate professor, Laboratory of Mechanics of Technological Processes

Russian Federation, 101-1, Prospect Vernadskogo, Moscow, 119526, Russian Federation.

Montaser Fekry

Moscow Institute of Physics and Technology

Email: montaser.fekry@yahoo.com
ORCID iD: 0000-0002-5825-931X
9, Institutskiy per., Dolgoprudny, Moscow Region, 141701, Russian Federation

References

Supplementary files