Расчет температурного поля пластины при электронно-лучевой сварке Текст научной статьи по специальности «Математика»

УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА Том 148, кн. 4 Физико-математические пауки 2006

УДК 519.6

РАСЧЕТ ТЕМПЕРАТУРНОГО ПОЛЯ ПЛАСТИНЫ ПРИ ЭЛЕКТРОННО-ЛУЧЕВОЙ СВАРКЕ

М.А. Игнатьева, Р.Ф. Кадыров, А.Б. Мазо

Аннотация

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

Введение

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

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

• построение единой во всей области неравномерной сетки, сильно сгущающейся в окрестности сварного шва [1];

•

в области высоких градиентов температуры и перестраивающихся при движении луча [2|:

держит траекторию источника [3. 4]:

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

Одним из наиболее распространенных в инженерной практике коммерческих пакетов, предназначенных для моделирования термомеханических процессов, является MSC.Marc-Mcrit.at [6]. включающий в себя и средства построения конечноэлементных сеток. В основу математической модели тепловых процессов с фазовыми превращениями металла при сварке в данном пакете положено стандартное

903 844 786 728 669 611 553 494 436 378 320

Рис. 1. Копечпоэлемептпая сетка и температурное поле при моделировании электроннолучевой сварки в пакете МБС.Магс

квазилинейное уравнение теплопроводности с объемно-распределенным нестационарным источником. Решение соответствующей краевой задачи осуществляется на основе ее конечно-элементной аппроксимации с последующим применением итерационного метода. Пример неравномерной сетки, сгущающейся в окрестности траектории движения нагревателя, а также температурное поле, построенное средствами МБС.Магс. показано на рис. 1.

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

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

1. Постановка задачи

Рассматривается процесс электронно-лучевой сварки пластины из алюминиевого сплава, которая имеет форму параллелепипеда со сторонами Ьх = 8 см,

Ьу = 2 см, Ьх = 5 см в направлении осей х, у и г соответственно. Теплопроводность к и теплоемкость с металла зависят от температуры Т, как показано на рис. 2; плотность р считается постоянной и равной 2710 кг/м3. Потерн энергии на плавление и выделение ее при кристаллизации (скрытая теплота фазовых переходов) учтены в функции с(Т) в виде резких пиков.

к, Вт/(м*К)

с, КДж/кг3

Рис. 2. Зависимости теплопроводности к и теплоемкости с от температуры Т

Электронный луч действует под углом а = 40° к нижней грани и движется с постоянной скоростью V = 2.5 см/с парадлельно оси OZ. Тепловая мощность нагревателя в каждый момент времени г распределена вдоль линии

г*(г) = vt, ж*(у) = Х0 + y/tg а, (1)

где жо = 0.5 (Ьх — Ьу/tgа) (см. рис. 3).

Распределение тепловой мощности, выделяемой электронным лучом, в каждый момент времени £ можно задать функцией д(у), где у — расстояние от верхней грани до луча. Определение погонной мощности д(у) для конкретного способа электронно-лучевой сварки представляет самостоятельную проблему [7]. однако в нашей задаче она считается заданной. Типичный вид этой функции представлен на рис. 4.

Математически удельная мощность Q(x,y, г,£) Ът/м3, выделяющаяся на линии движущегося луча, может быть записана с помощью произведения 6-функций Дирака:

Q(x,y, г, г) = я(у)6[ж — ж*(у)]6[г — г*(г)],

где координаты ж* и г* заданы формулами (1). При проведении расчетов исполь-

6

( — совГ^У |£|<г,

<*(£) = <4г \2г)’ ’

[0, |£|> г.

q, МВт/м

Рис. 4. Погонная мощность источника д(у) в зависимости от расстояния до верхней грани

Здесь г - радиус регуляризации (в расчетах принималось г = 5 мм, что соответствует реальному радиусу теплового воздействия электронного луча). Окончательная формула для мощности источника в произвольной точке области х, у, г в момент времени Ь принимает вид

(2)

0,

,

где £ = х — х*(у), п = г — г*(Ь).

Температурное поле Т(х,у,г,Ь) в пластине, представляющей из себя параллелепипед О = (0, Ьх) х (0, Ьу) х (0, Ьг) с границей Г, описывается начально-краевой задачей

3Н (Т)

дЬ

— ё1у к(Т)grad Т = Q(x, у, г, Ь), (х, у, г) € О, Ь > 0,

дТ

к{Т)— = 0, (х, у, г) € Г, # > 0,

Т(х, у, г, 0) = Т°(х, у, г), (х,у,г) € О.

(3)

Здесь п — внешня я к Г нормаль, движущийся объемно-распределенный источник ^ ^^^улой (2), а Н - удельная энтальпия, Дж/м3, определяемая как

Т

Н (Т) = рI е(0<%.

То

Математическую формулировку (3) назовем нестационарной задачей. При этом среда (металл свариваемого изделия) покоится, а источник (электронный луч) движется. В ряде теоретических работ по моделированию сварочных процессов данная постановка называется задачей в переменных Лагранжа [8, 9]. Именно такая математическая модель реализована в пакете МБС.Магс.

Альтернативный подход к математическому описанию изучаемого процесса состоит во введении системы координат, связанной с движущимся нагревателем (подход Эйлера, [9]). В этом случае источник тепла становится неподвижным, а

г

холодный металл набегает на него со скоростью сварки v. При v = const в дифференциальном уравнении теплопроводности исчезает производная по времени, но появляется конвективный член: дН

V— - div к(Т)grad Т = Q(x, у, z, t), (ж, у, zj G П,

Т(x,y,z) = To, (x,y,z) G Го (4)

дТ

k(T)— = 0, (ж, у, z) G Гдг U Ть,

где rN = {(x,y,z) G Г : z G (0, Lz)}, Го = {(x,y,z) G Г : z = 0}, Гь = {(x,y,z) G G Г : z = Lz}. В постановке (4), которую будем называть стационарной задачей, время t в правой части уравнения является параметром и служит лишь для определения положения источника по формулам (1).

2. Численное решение задачи

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

Проведем неявную полудискретизацию уравнения, заменив частную производную дН/дЬ разностным отношением (Н — Н)/т, где Н = Н(Ь — т), т — шаг по времени. В результате на каждом временном слое получим нелинейное уравнение А Т = Е, где

АТ = Н(Т) — тdivк(ТТ, Е = Н + тQ.

После пространственной дискретизации этого уравнения получим систему нелинейных алгебраических уравнений относительно узловых значений температуры Т^, г = 1,...,М. Эту систему будем решать итерационным методом релаксационного типа, в котором для каждого г требуется решить скалярное нелинейное уравнение вида А(Т) = ¥1. Его единственное решение всегда существует, поскольку А(Т) - монотонно возрастающая функция.

Численное решение этого уравнения может быть получено быстро, если вместо А( Т)

ризуем пространственный оператор div к(Т)grad Т, введя функцию Кирхгофа

Т

К(Т) = ! к(£Щ. (5)

То

Тогда определяющее уравнение примет вид

Н(Т) — Н(Т) _ = ^

т

где Д - оператор Лапласа. Для уравнения (6) функция А(Т) имеет вид А(Т) = = аН(Т) + вК(Т), поэтому для ее кусочно-линейного представления достаточно в том же виде представить функции Н(Т) и К (Т).

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

Для стационарной задачи (4) также проводится линеаризация (5), в результате получается уравнение конвективной теплопроводности

- АЩТ) = д, (7)

которое после пространственной дискретизации также решается методом релакса-

Применение итерационного метода предполагает задание начального приближения. В нестационарной задаче (6) выбор хорошего начального приближения очевиден: это значение температуры на предыдущем временном слое. Однако, для того чтобы задать начальное приближение, близкое к решению стационарной задачи (7). требуются дополнительные усилия. Рассмотрим уравнение (7) с постоянными коэффициентами к и с, определенными как интегральные средние от к и с то температуре. В этом случае Н = с рТ, К = кТ, и уравнение становится

линейным: УдТ/дг — ДТ = /, где V = 'ирс/к, / — Q/k. Решение этого уравнения

будем искать в виде

Т (х) = и(х)ввг, в = У/2. (8)

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

Уи — Ди = / = в-вг/, (х, у, г) £ П;

ди

-----/Зи = 0, (х,у, г еГь;

дп (9)

ди

7^=0, [х,у, г) £ Гдг;

и(х, у, г) = в-вгТо, (х, у, г) £ Го.

Оператор задачи (9) линеен, симметричен и положительно определен, поэтому для ее численного решения применимы известные эффективные итерационные методы.

и

Т

в качестве начального приближения в методе релаксации для задачи (7).

3. Триангуляция области и сеточные аппроксимации

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

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

Для создания исходной сетки построим вначале в сечении ХУ плоскую сетку четырехугольиых элементов, топологически эквивалентную ортогональной, со сгущением узлов к линии действия источника у = tg а(х — Ьх/2) — Ьу/2. Пример такой сетки показан на рис. 5. Построение трехмерной сетки на первом этапе за-

Рис. 5. Сетка в сечении ХУ

Рис. 6. Композиционная сетка в сечении ХУ (справа увеличенный фрагмент)

г

отрезка [0,Ьг].

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

ХУ

тырехуголышка снова строится аналогичная сетка с тем же законом сгущения и с одинаковым количеством «мелких» элементов внутри каждого «крупного».

ХУ

ставлен на рис. 6. Построение мелкой сетки в выделенной подобласти завершается

г

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

45

40

35

30

25

20

15

10

5

0

г, мм 49

48

47

46

45

44

43

10 20 30 40 50 60 70 х, мм

13 14 15 16 17 18 19 20 21 х, мм

Рис. 7. Композиционная сетка в сечении XZ (справа увеличеннвш фрагмент)

мм

Для стационарной задачи точность численного решения в значительной степени определяется способом аппроксимации конвективного члена. Известно, что направленная разность «против потока» имеет первый порядок аппроксимации и вносит в сеточную схему регулярное возмущение в виде дополнительной численной вязкости (в нашем случае - теплопроводности). Аппроксимация первой производной центральной разностью избавлена от этого недостатка, однако приводит к плохо обусловленной матрице системы сеточных уравнений. Одним из способов повышения порядка аппроксимации конвективного слагаемого в задаче (4) при сохранении диагонального преобладания матрицы соответствующей системы уравнений может быть его представление в виде

г>^«г>я5+г>(у-е)яа5, (10)

где е £ (0, Нг/2] - варьируемый параметр.

Данная аппроксимация была применена для стационарной модели со значениями е = 10-3Нг, Нг/2. Проводилось сравнение полученного решения с решением нестационарной задачи в соответствующий момент времени. Без регуляризации (е = Нг/2) наблюдалось значительное (порядка 150 К) расхождение решений, тогда как при е = 10-3 Нг эти решения практически совпал и (кривые 1 и 2 па рис. 8).

4. Численные эксперименты

Описанные алгоритмы были численно реализованы на языке С--. Размерности

сеток были выбраны следующими. Количество узлов крупной сетки в О составило 41 х 21 х 41 то направлениям х, у и г соответственно, а мелкой сетки в подобласти

— 33 х 61 х 121. Шаг то времени т для нестационарной модели равнялся 0.01 с.

Система нелинейных уравнений во всех случаях решалась методом верхней релаксации с итерационным параметром ш = 1.6. Численное решение стационарной задачи проводилось при двух разных положениях источника: г0 = Ьг/4, Ьг/2. На рис. 9 изображены изотермы на верхней грани и в плоскости действия источника, построенные по результатам решения нестационарной задачи. Картина температурного поля качественно согласуется с решением нестационарной задачи (3), полученным в пакете МБС.Магс на неравномерной сетке с 13849 узами и 12160 элементами. Наибольшее отличие наблюдалось в окрестности действия источника 3

т, К т, К

Рис. 8. Температурные кривые вдоль траектории движения источника при при ,г* = Ьх/4 (слева) и при = Ь/2 (справа)

г, мм 45 40 35 30 25 20 15 10 5 0

0 10 20 30 40 50 60

у, мм

15

10

300 400 500 600 700 800 900

300 400 500 600 700 800 900

Рис. 9. Температурное поле па верхней грани (слева) и в плоскости действия источника (справа) при ,г* = Ь/2

Очевидно, что нестационарная модель (3) более точно описывает нагрев конечного тела движущимся источником, поскольку стационарная постановка (4) справедлива, строго говоря, лишь в бесконечной области. Количественное влияние границ на решение стационарной задачи иллюстрирует рис. 10. Нетрудно видеть, что при расположении нагревателя в средней части траектории оба решения близки. а когда центр источника приближается к границам пластины на расстояние порядка г (см. (2)), погрешность стационарной модели существенно возрастает.

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

• “НСМ” — нестационарная модель (3) с использованием композитной сетки и метода верхней релаксации:

•

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

(9). которое находится методом сопряженных градиентов:

Результаты численных экспериментов сведены в табл. 1. где указаны время расчета н среднеквадратичная норма отклонения от эталонного решения. В качестве эталонного принималось решение “НСМ” в соответствующий момент времени.

|Т, - Т 21, К

Рис. 10. Максимальная разность между решениями стационарной и нестационарной задач в зависимости от положения источника

Табл. 1

Характеристики вычислительной эффективности различных методов расчета температурного поля

г* МвС.Маге НСМ СМ СМНП

ьг £. мил 10.05 1.66 1.02 0.99

4 ||Т-Т*|| 30 0 0.098 0.098

ьг £. мил 19.97 3.288 1.172 1.141

2 ||Т-Т*|| 51.3 0 0.027 0.027

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

Выводы

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

2. Показано, что стационарная задача с достаточной точностью описывает динамику температурного поля, когда электронный луч расположен на расстоянии более 0.5 см от границ пластины (рис. 10).

3. Сравнительный анализ эффективности различных методов показал, что быстрее других работает алгоритм решения стационарной задачи конвективной теплопроводности со специальным выбором начального приближения, однако он дает значительные ошибки решения в начале и конце сварки. Незначительно уступает ему по скорости метод решения задачи в нестационарной постановке с движущимся нагревателем. Оба предложенных алгоритма работают значительно быстрее. чем алгоритм, реализованный в МБС.Магс: расчет на композиционной сетке, содержащей в 20 раз больше узлов, чем сетка в МБС.Магс. требует приблизительно в 6 раз меньше времени (табл. 1).

Работа выполнена при финансовой поддержке гранта INTAS-Airbus (проект .V 04-80-6951).

Summary

М.A. Ignatieva, R.F. Kadyrov, А.В. Mazo. Calculation of the temperature field of a plate when elect.ron-beam welding.

Numerical methods for solving of temperature field dynamics problem for elect.ronic-beam welding are proposed. Stationary and lion-stationary formulations of the problem are considered. In a FEM approximation, composite mesh with a local refinement around the weld seam is used. Numerical efficiency of constructed methods is analysed as well as the comparison of these methods with MSC.Marc solver is executed.

Литература

1. Bangerth W., Rannacher R. Adaptive finite element methods for differential equations.

Basel: Birkliauser Verlag, 2003. 207 p.

2. Гильманов A.H. Методы адаптивных сеток в задачах газовой динамики. М.: Физ-

матлит. 2000. 248 с.

3. Quarteroni A., Periaux J., Kuznetsov Yu.A., Widlund О.В., Eds. Domain decomposition

methods in science and engineering. The Sixth Int. Conf. on Domain Decomposition. V. 157 of Contemporary Mathematics. Providence. RI: AMS, 1994. 485 p.

4. Агоитов В.И. Методы разделения области в задачах математической физики // Вычислительные процессы и системы. М.: Наука, 1991. Вып. 8. С. 4 51.

5. Баранов П.А., Исаев С.А., Кудрявцев Н.А., Харченко В.Б. Расчет колебаний цилиндрического маятника в наполненной вязкой жидкостью полости с использованием скользящих мпогоблочпых сеток // ИФЖ. 2003. Т. 76, Л'! 5. С. 61 70.

6. MS С. Software Corp. MSC.Marc: Theory and User Information // MSC.Software Corp.: Ver. 2005, printed in U.S.A.

7. Судник B.A., Карпухин E.B., Padau Д., Хекелер Г. Метод эквивалентного источника теплоты: поверхность солидуса движущейся сварочной ваппы для МКЭ-расчета трехмерного температурного поля // Компыот. техп. в соед. материалов. Тр. 2-й Всерос. пауч.-техп. копф. Тула: ТулГУ, 1999. С. 49 63.

8. Shanghvi J.Y., Michaleris P. Tliermo-elasto-plastic finite element analysis of quasi-state processes in Eulerian reference frames // Int. J. Numer. Met.li. Eng. 2002. V. 53.

P. 1533 1566.

9. Zhang Lu, Michaleris P. Investigation of Lagrangian and Eulerian finite element methods for modelling the laser forming process // Finite elements in analysis and design. 2004. V. 40. P. 383 405.

10. Ewing R.E., Lazarov R.D. Local refinement techniques in the finite element and finite difference methods // Proc. Int. Conf. of Numerical Methods and Applications. Sofia: Publishing House of the Bulgarian Academy of Sciences, 1989. P. 148 159.

11. Kadupoe Р.Ф., Мазо А.Б. Расчет тепловых полей при сварке пластип движущимся источником // Тр. XI Всеросс. шк.-семинара «Современные проблемы мат. моделирования». Ростов п/Д: Изд-во Ростовск. ун-та, 2005. С. 161-170.

Поступила в редакцию 04.10.06

Игнатьева Марина Александровна кандидат физико-математических паук, старший научный сотрудник НИИ математики и механики им. Н.Г. Чеботарева Казанского государственного университета.

Е-шаП: Marina.IgnatievaQksu.ru

Кадыров Рафаэль Фаридович научный сотрудник НИИ математики и механики им. Н.Г. Чеботарева Казанского государственного университета.

Е-шаП: гкайугт Qksu.ru

Мазо Александр Венцианович доктор физико-математических паук, профессор кафедры аэрогидромехапики Казанского государственного университета.

Е-шаП: атаги Qk.su. ги