CC BY
49
Вычислительные технологии
Том 11, № 4, 2006
КОНЕЧНО-РАЗНОСТНЫЙ АЛГОРИТМ ДЛЯ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ ПРОЦЕССОВ ЛАЗЕРНОЙ СВАРКИ МЕТАЛЛИЧЕСКИХ
ПЛАСТИН*
В. П. Шапеев, А. Н. Черепанов Институт теоретической и прикладной механики СО РАН,
Новосибирск, Россия e-mail: [email protected], [email protected]
A finite-difference algorithm is formulated on the basis of thermophysical model for numerical simulation of laser butt welding of metal plates. Numerical simulation of welding of aluminium alloy plates is carried out. All physical constant needed for the calculations are taken from the known sources. In the problem of forming weldpool without steam channel, influence of heat exchange increase in the weldpool on the weldpool size is investigated. It is shown that in absence of steam channel essential increase of effective heat conductivity inside the pool leads to relatively small increase in its size and thickness of the welding plates.
Введение
Создание эффективной технологии лазерной сварки — одна из актуальных задач современных авиастроительной, космической и других отраслей промышленности. Лазерная сварка имеет ряд преимуществ перед другими широко используемыми способами соединения деталей. Однако в настоящее время она нуждается в разностороннем совершенствовании, в частности в повышении стабильности и управляемости процесса. Во многих странах с развитой промышленностью активно ведутся исследования в этом направлении.
Сложность физических процессов, протекающих при высокой температуре в малых по размерам зонах воздействия лазерного луча на материал, затрудняет экспериментальные исследования лазерной сварки. Математическое моделирование может дополнить и частично заменить физические эксперименты. В этом направлении ведется активный поиск, и отдельные параметры процесса поддаются расчетам. Но на данный момент нет завершенной модели, которая описывала бы во всех взаимосвязях основные явления, определяющие этот сложный физический процесс [1-14]. В частности, при достаточно высокой мощности лазерного излучения в зоне воздействия луча на свариваемый материал образуется паровой канал. При этом материал свариваемых деталей (чаще всего металл)
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 06-01-00080-а).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.
находится в разных агрегатных состояниях. Собственно процессы в зоне парового канала, скорость и характер их протекания, температура и другие параметры определяют устойчивость процесса сварки и качество сварного шва. Здесь необходимо учитывать сложные взаимодействия ванны как с жидким металлом, в котором образовался паровой канал, так и с остывающим металлом.
В настоящей работе излагается численное моделирование процесса лазерной сварки встык металлических пластин на основе описанной в [14] теплофизической модели. В отличие от [14] здесь приводится численный алгоритм решения задачи, а сведения по физической модели даны в минимальном объеме, необходимом для понимания численного алгоритма и полученных результатов.
1. Постановка задачи
Пластины в форме прямоугольных параллепипедов одинаковой толщины состыкованы между собой узкими боковыми гранями. Поверхности пластин для предохранения во время сварки металла от нежелательных окислительных процессов обдуваются инертным газом. Ось луча лазера лежит в плоскости стыка пластин и направлена перпендикулярно к их поверхности. На рис. 1 схематично показана область сварки.
Начало декартовой системы координат находится на пересечении оси луча лазера с верхней поверхностью пластин. Ось г направлена вдоль оси луча, ось х — в плоскости стыка пластин и направлена в сторону, противоположную направлению перемещения луча лазера, а ось у перпендикулярна к стыку. Луч лазера движется со скоростью (сварки) V против направления оси х. В выбранной системе координат, движущейся вместе с лучом, пластины перемещаются со скоростью V, а луч неподвижен. Для наглядности на рисунке приведен разрез свариваемой пластины с удаленной частью в первом квадранте плоскости х, у. Паровой канал 2, формирующийся в окрестности оси луча 1, граничит с подобластью 3, занятой жидким металлом. Последняя граничит с подобластью 4, занятой двухфаз-
Рис. 1. Схема сварочной ванны (сечение у = 0): 1 — луч лазера; 2 — парогазовый канал; 3 — жидкая фаза; 4 — двухфазная зона; 5 — твердая фаза.
1
2
ным состоянием вещества, которая в соответствии с физикой процесса имеет небольшую толщину и контактирует с периферией зоны сварки 1, занятой твердым металлом. Топология названных свободных криволинейных границ (поверхности парового канала и границ области с двухфазным состоянием металла) зависит от мощности лазера.
В математической модели процессы теплопереноса при сварке описываются квазистационарным уравнением теплопроводности, а теплообмен с окружающей средой и между подобластями с металлом в различных фазах — различными нелинейными граничными условиями теплового баланса. Коэффициент температуропроводности — константа, имеющая в разных подобластях расчетной области значения, равные некоторым средним его значениям для конкретной фазы металла в соответствующей подобласти.
Для упрощения этой сложной задачи трехмерное квазистационарное уравнение теплопроводности осредним по координате y, в результате чего получим
дТ _ Л f д2Т i д2Т \ \г где эффективная теплоемкость сег вычисляется по формуле
= Ч + - ж(T - (1)
ClPl, T < Te
e j
Cei = < 1 + C2, Te < T < Tío, (2)
СзРз, Ti o < T.
Здесь индекс i = 1, 2, 3 обозначает параметры твердой, жидко-твердой и жидкой фаз соответственно; c¿, p¿, A¿ — теплоемкость, плотность и теплопроводность i-й фазы соответственно; к — удельная теплота плавления; Ti0, Te — температуры начала и конца затвердевания; it = 2^/a¡r — длина распространения тепловой волны за время т = 2tf/v, tf — радиус луча в фокальной плоскости; a¿ — температуропроводность i-й среды; Ta — характерная температура пластины вне области осреднения. Как вычисляется сечение (доли) жидкой фазы fi (T) в двухфазной зоне, описано в [14].
Граничные условия задачи получаются из условий теплового баланса. На поверхности z = 0 в областях твердого, затвердевшего, а также жидкого сплава и двухфазного состояния условие имеет вид
дТ
dz
z=0
= (ak + )(T |z=o - Tgi ) , г = 1, 2, 3. (3)
Здесь ак — конвективный коэффициент теплоотдачи [9], который обусловлен обдувом зоны сварки инертным газом и зависит от характера его течения (ламинарное или турбулентное); Тд — температура газа; ап - радиационный коэффициент теплоотдачи, определяемый соотношением
аГг = егаа (Т 12=0 + Тд2) (Т |*=, + Тд), (4)
где а0 — приведенная степень черноты и константа Стефана — Больцмана соответственно.
Условие, аналогичное (3), имеет место на нижней поверхности свариваемых пластин:
dT
dz
= (ak2 + an ) (T |z=h - Tg2 )
z=h
где ап = £-¿00 (Т || + Тд22) (Т + Тд2), Тд2 — температура среды (монтажного стола), контактирующей с нижней поверхностью пластин; — коэффициент теплоотдачи от нижней поверхности к монтажному столу. При контактном охлаждении через зазор, заполненный газом, ак2 = Ад/8д (Лд, 8д — теплопроводность газа и толщина газового зазора между нижней поверхностью пластины и поверхностью монтажного стола). В случае, когда и нижняя поверхность пластины обдувается газом, вместо следует брать зависящий от характера течения коэффициент ак.
Условие теплового баланса в области действия лазерного излучения (х < г^) на поверхности парового канала Zc(x) имеет вид
—А3УТ ■и = дп — Ьт + 5дс. (5)
Здесь т — массовая скорость испарения вещества с единицы поверхности, связанная с избыточным давлением паров Р(г), необходимым для удержания стенок канала от схло-пывания, соотношением т = л/Р(г)рь, где р.и — плотность пара; Ь — удельная теплота испарения сплава; д — поглощенный поток излучения с учетом переотражения; 8 = 0 на передней стенке канала, 8 =1 — на задней. Величина qc — поток энтальпии, переносимый из области передней стенки канала к задней, определяемый соотношением qc = о3р^г (Т1—Т), где Т1 — среднее значение температуры на передней стенке канала; vr — средняя скорость движения расплава, обтекающего стенки канала (в процессе перемещения канала вдоль оси х), значение которой оценивалось из уравнений расхода и неразрывности [6, 10].
Для интенсивности излучения лазера принимаем нормальный закон распределения Гаусса, и оно зависит от переменного по г радиуса луча лазера
7 . 2\ 1/2
2 ¡г — Zp
Гг = \ГР + -Ао
\ ПГр
Здесь ZF — координата по г фокальной плоскости, г^ — радиус лазерного луча в фокальной плоскости, А0 — длина волны излучения лазера (в приведенных ниже расчетах А0 = 10.6 мкм). Считаем, что остаток мощности после первого отражения (Ш1 = (1 — А) Ш, где А — коэффициент поглощения) равномерно рассеивается по поверхности канала Бс, состоящей из боковой поверхности канала Бъ, входного и выходного отверстий.
Мощность, поглощаемую стенками канала, вычислим приближенно через эффективный коэффициент поглощения Аец = в А + Ае(1 — А) по формуле ШаЪз = Ае^Ш (в =1 на пересечении Бс с лучом лазера и в = 0 вне луча). Эквивалентный коэффициент поглощения для многократно отражаемых излучений лазера определяем по формуле [3]
Ае = А£
п= 1
(1 — А) Бъ
Бс
(6)
Тогда в (5) плотность мощности излучения лазера, поглощенную элементом поверхности парогазового канала г = Zc(x,y), найдем как
q(x, y, Zc(x, у)) = ехр (-2г2/г2) . (7)
пг2
п
Далее проведем численное моделирование сварки металлических пластин на основе конечно-разностного метода. Поскольку положение поверхности канала и границ между
указанными подобластями зависит от неизвестного решения и краевые условия нелинейны, вычислительный алгоритм будет итерационным с вложенными циклами. Для удобства программной реализации вычислительного алгоритма и его изложения все множество итераций во внешнем цикле разбито на этапы. Расчетная область О на первом этапе вычислений (без парового канала и валика из затвердевшего металла) — прямоугольник в плоскости стыка пластин. Его стороны г = 0, £ = к, х = -1\, х = /2, где 1\ и /2 — расстояния соответственно от левой и правой границ прямоугольника до оси луча лазера. Выбор величин 1\ и 12 зависит от мощности лазера и скорости сварки. При этом основное требование такое, чтобы на левой ("невозмущенной") границе с хорошим приближением можно было значение температуры полагать равным температуре внешней среды ("на бесконечности"), т. е. в качестве краевого условия на этой границе принимаем
В данных расчетах при выборе краевого условия на правом краю (x = l2) прежде всего учитывалось поведение решения данной задачи. Заметим, что для случая, когда толщина пластины (ширина расчетной области) много меньше ее длины (длины расчетной области), по мере удаления вправо от области жидкой ванны с паровым каналом температура и ее производная по x выравниваются по толщине пластины (сечению x = const). Пусть длина этой области такова, что на правом краю находится достаточно протяженная подобласть с твердым металлом. В этом случае постановка условия Дирихле с температурой T'х или некоторой другой, меньшей температуры плавления, мало влияет на распределение температуры в интересной для нас зоне парового канала и ванны. Значение температуры на правой границе длинной области существенно влияет только в подобласти, прилегающей к этой границе. Для выбора приемлемого значения температуры на правой границе решалась прикидочная модельная задача со всеми условиями на границах и конкретными значениями физических параметров, с которыми решалась основная задача. Но при этом в отличие от варианта основного расчета область бралась прямоугольной (без постановки парового канала) и значительно более протяженной вправо. На правой границе ставилось условие Дирихле со значением температуры Полученное распределение температуры позволило установить некоторое приближенное распределение температуры (и ее градиента) по x в промежуточной зоне, находящейся на некотором удалении вправо от ванны и влево от удаленной правой границы. Этот расчет брался за основу при выборе значения температуры в условии Дирихле на правом краю. О точности выбора можно судить по косвенному признаку: насколько соответствуют друг другу распределение температуры вблизи правого края и градиенты температуры при (x = l2) в расчетах по упрощенной и более полной моделям.
Разумеется, при постановке условия на правом краю подходы, основанные на учете потоков тепла через все внешние границы расчетной области, и требования выполнения полного баланса тепла имеют шанс быть более точными. Но они, во-первых, трудно осуществимы, во-вторых, сложно доказать сходимость некоторых итерационных процессов, которые при этом возникают.
(8)
2. Аппроксимация уравнения
Область решения задачи покрыта прямоугольной разностной сеткой О^, с узлами (гкх, , г = 0,1,..., П\, ] = 0,1,..., п2, где Н\ и к2 — шаги сетки. Для численного решения задачи
использованы варианты схемы метода установления. Одна из них для уравнения (1) имеет вид
тп+1 _ ТП ТП+1 _ ТП+1
г,3_Ч? + г+1 _
Т 2^ _
/ Tn+1 _ отn+1 I Tn+1 Tn _ OTn i T™ i
ai-Г2- + -Г2- - 72(1i,i — 1a)■
/¿1 /^2 /
Кроме нее была испытана также схема условной аппроксимации. Отметим, что сходимость численного решения наблюдалась по сходимости величины max |Tin7+1 — Щ | в расчетной
области. При использовании второй схемы дополнительно наблюдалась величина невязки от подстановки разностного решения в разностное уравнение теплопроводности и в разностные условия на границах. При этом, естественно, в них удалялись фиктивные производные по времени, а производные по пространственным переменным и неоднородные члены все выписывались на последнем временном слое. Итерации прекращались, как только эти величины становились меньше некоторых заданных величин. Численное решение с точностью погрешности аппроксимации по обеим схемам получилось одно и то же. Поэтому можно говорить, что вторая схема использовалась для контроля численного результата, полученного по первой схеме.
Схема (9), неявная по x и явная по z, реализуется итерационно прогонками по индексу i (т.е. вдоль координатных линий x = const). Таким образом, для каждого значения j прогонкой решаются системы уравнений с трехдиагональными матрицами.
3. Аппроксимация краевых условий и связанные с ними некоторые детали вычислительного алгоритма
В данной работе краевые условия на внешней и внутренних границах (фронтах между фазами материала изделия) аппроксимированы с первым порядком. Для записи краевых условий в конкретной точке необходимо иметь численные значения компонент нормали, кривизны границы. Поскольку все внутренние границы в рассматриваемой задаче — это определенные изотермы, при необходимости по вычисленным сеточным значениям температуры интерполяцией находятся соответствующие границы. Граница между жидким металлом и двухфазной зоной — изотерма г _ Zl(ж), которая определяется по температуре ликвидуса сплава Тю. Иначе говоря, Тю — это значение температуры, отделяющее в процессе остывания жидкое состояние металла от его двухфазного состояния. Граница между двухфазной зоной и твердым металлом г _ Ze(ж) — изотерма, которая определяется по температуре эвтектики Те — температуре в конце затвердевания жидкой фазы в процессе остывания. В нашем рассмотрении указанные изотермы на границе расчетной области (г _ 0) имеют по одному концу в области положительных (со стороны задней стенки парового канала) и отрицательных (со стороны передней стенки парового канала) значений ж. В зависимости от толщины пластин и величины парового канала каждая из этих изотерм может состоять из одного куска либо из двух разъединенных кусков. Во втором случае они имеют концы на границе области г _ Л. При этом во всех расчетах, проведенных в данной работе со входными данными, взятыми из практики сварки, во втором случае на рассматриваемых кривых зависимость ж _ ж(г) была однозначной. А в
первом случае всегда имелась точка экстремума по z, которая отделяла на каждой кривой две ветки с однозначной зависимостью x = x(z). Это справедливо и по отношению к линии z = Zc(x), моделирующей поверхность парового канала. Указанное важное свойство использовалось при дискретной аппроксимации этих границ и условий на них.
Процесс построения изотерм автоматизирован и проводится по однозначному алгоритму, подробное описание которого здесь излишне. Назовем только отдельные его детали. На разностной сетке путем сравнения значений температуры на координатных линиях (x = const и y = const) определяются две соседние точки, между которыми находится конкретная изотерма. Приближенно определяется ее угол наклона к оси x. Затем она разбивается на отдельные дуги, синус угла наклона которых к оси x меньше \/2/2 (пологие участки), и дуги с углом наклона ближе к вертикальному (крутые участки). Положение точек на пологих участках изотерм определяется интерполяцией по значениям температуры в соседних узлах на сеточной линии y = const, между которыми находится изотерма. Положение точек на крутых участках находится аналогичной интерполяцией на сеточных линиях x = const.
На изотерме ликвидуса z = Z (x) требуются непрерывность температуры и условие баланса тепла
Л —I = л —I (10)
dn lz—Zi dn
На изотерме эвтектики г = Д=(х) требуются непрерывность температуры и условие баланса тепла
Je
Л2 Tl—z+ = Лз Z- + PeSt . =. (11)
dn lz~Ze dn12-Ze I f dz 2
dx
В (10) и (11) индексы " +" и " —" означают, что соответствующие производные от температуры взяты на кривой как предельные с разных ее сторон. Последнее слагаемое в (11) связано с теплотой кристаллизации металла. Коэффициентом при нем является произведение чисел Пекле и Стефана.
Для того чтобы аппроксимировать соотношения (10) и (11), на сеточном решении задачи найденные линии г = ^(х) и г = Д=(х) аппроксимируются ломаными линиями, проходящими через узлы сетки. Для этого из каждой пары узлов сетки, между которыми находится рассматриваемая линия, выбираются ближайшие к ней. Выбранные узлы можно последовательно и непрерывно соединить между собой отрезками из соответствующих координатных линий или диагоналей ячеек сетки так, что на полученной ломаной не будет прямых углов. В представленной работе при аппроксимации величин дТ/дп в условиях на внутренних (и внешних) границах производные дТ/дх, дТ/дг аппроксимировались с первым порядком с помощью значений температуры в узлах сетки с соответствующей стороны рассматриваемой линии. Заметим, что концы этих линий находятся на внешней расчетной границе области и на них выписываются условия на внешней границе.
Для ускорения итерационного процесса во все разностные граничные условия искусственно добавлялась фиктивная разностная производная по времени от температуры. Так же как и в уравнении (9), знак при ней выбирался таким, чтобы увеличивался модуль элемента на диагонали в соответствующей строке матрицы системы линейных алгебраических уравнений разностной задачи.
Процесс численного построения канала состоял из нескольких этапов. На первом этапе задача численно решалась в прямоугольной области без парового канала. Итерационный
процесс начинался, когда температура во всей области была равна температуре окружающей среды. После того как максимум невязки по всем узлам сетки от подстановки численного решения в разностные уравнения становился меньше заданной величины, согласно описанной выше процедуре на сетке находились изотермы кипения металла, ликвидуса и эвтектики. Тем самым находилось некоторое грубое приближение решения и границ, используемое в качестве начального для решения на последующих этапах его уточнения.
В данной работе выбран способ итерационного процесса, в котором решение задачи с нелинейными краевыми условиями третьего рода на отдельной итерации по всей расчетной области сводилось к решению системы линейных алгебраических уравнений. Он состоит из нескольких известных приемов. Так, проводится линеаризация нелинейных уравнений. В краевом условии (3) нелинейный коэффициент ar1 аппроксимируется на решении, вычисленном на предыдущей итерации.
На втором этапе полученная расчетная область разбивается на подобласти. Они отличаются между собой либо коэффициентом теплопроводности, либо сочетанием краевых условий на двух концах каждой координатной линии x = const, упирающейся концами в границы конкретной подобласти. Сеточные линии x = const внутри этих подобластей не пересекаются их границами. Число подобластей зависит от наличия парового канала и от того, как расположились линии ликвидуса и эвтектики.
В процессе проведения одной текущей итерации во всей расчетной области (глобальной итерации) система (9) прогонками вдоль оси x последовательно решается во всех подобластях. Внутри подобласти последовательно перебираются все возможные в ней значения индекса j. Значения решения на границах всех подобластей берутся с предыдущей итерации. После завершения итерации уточняются значения решения во всех узлах на границе подобластей. В каждом узле на границе краевое условие разрешается относительно значения температуры в выбранном узле. Тем самым находится его уточненное значение через значения температуры в соседних узлах, полученные на данной итерации. Оно и используется на следующей итерации в качестве значения решения на границе.
Совокупность описанных приемов позволила получить достаточно простой способ реализации численной модели.
4. Моделирование парового канала с использованием условия теплового баланса
Построение парового канала и моделирование физических процессов на его границе — это самая сложная часть решения рассматриваемой задачи. Отметим, что поведение решения на границах канала играет ведущую роль в итерационном процессе. Решение устанавливается прежде всего в окрестности парового канала, а затем во всей расчетной области. От размера и формы парового канала зависят форма и положение изотерм в расчетной области.
Здесь уместно отметить, что на первом этапе решения задачи коэффициент теплопроводности во всей области брался соответствующим жидкому металлу. Нетрудно понять, что при решении задачи в рамках рассматриваемой физической модели максимальное значение координаты z, полученное на изотерме кипения металла z = Zsat(x), служит некоторой оценкой снизу глубины парового канала. Если она оказывается больше толщины пластин, то ясно, что при данной мощности лазера паровой канал пробивает их насквозь. В этом случае, видимо, лазер данной мощности можно использовать для резки
пластины заданной толщины из данного металла.
Естественно, после первого этапа есть оценка снизу ширины парового канала. Ее определяет положение концов изотермы кипения на поверхности (г = 0). Они могут служить в качестве ориентира при выборе приближенной формы канала на втором этапе.
Если форма границ расчетной области известна, то задача со сформулированными граничными условиями имеет решение и строится численно. Однако форма парового канала неизвестна. Вопрос строгого выяснения факта существования решения задачи, которая здесь возникает, не является предметом данной работы.
Форма парового канала зависит прежде всего от характера процесса испарения жидкого металла с его поверхности, давления пара, силы поверхностного натяжения и гидростатической силы. Величина последней сравнительно мала для тонких пластин. Видимо, ее учет мало влияет на форму канала и в итоге на температурное поле. Давление пара складывается из давления насыщенного пара и силы отдачи от истечения пара из канала. Сила поверхностного натяжения (пропорциональная кривизне поверхности) стремится сделать плоской искривленную поверхность. В целом она стремится вернуть поверхность канала в плоскость поверхности пластин (г = 0) (кривизна которой равна нулю). Таким образом, сила поверхностного натяжения уравновешена давлением пара в канале. Давление пара зависит от процесса испарения жидкого металла. Более того, если бы отсутствовали давление пара и сила поверхностного натяжения, то форма канала полностью определялась бы только процессом испарения, которое учитывается здесь в граничном условии (5).
На втором этапе при построении границ парового канала используется поле температур, полученное на первом этапе. Форма передней границы канала (расположена перед движущимся лучом лазера) подбирается так, чтобы на нем прежде всего выполнялось краевое условие (5). Если при сварке кроме лазера нет других источников энергии, то наибольшая часть потока тепла извне внутрь области проходит через небольшую окрестность оси луча лазера на передней стенке канала. При этом температура, которая достигается здесь, зависит от угла наклона стенки. В регулировании величины угла наклона стенки канала к оси х определяющую роль играет следующий физический процесс. Поскольку величина главной части поглощенной энергии стенкой пропорциональна косинусу угла ее наклона (к оси х), то (при мощности лазера выше критической) при уменьшении угла наклона начинает испаряться больше металла, и канал углубляется, что ведет к увеличению угла наклона стенки. Однако это увеличение угла ограничено сверху. Во-первых, при недостаточной толщине пластины произойдет сквозной пробой пластины лазером, что нежелательно. Во-вторых, увеличение угла больше некоторого, соответствующего заданным параметрам задачи (мощности лазера, материала пластин и пр.), уменьшает количество поглощенной энергии через единицу площадки стенки, температура на стенке канала падает, а глубина канала уменьшается. Действие указанного механизма описывают первое и второе слагаемые в правой части уравнения (5). Третье слагаемое несоизмеримо мало по сравнению с ними. Особо отметим, что распределение энергии по радиусу в сечении луча представляет собой ярко выраженный пик с максимумом на оси луча. Поэтому наибольшая часть передней стенки канала в направлении оси г находится в окрестности оси луча. Более того, при скоростях сварки V, близких к оптимальным, приближенно полагаем, что передняя стенка на верхней поверхности пластин (г = 0) начинается в точке (х = — гр, г = 0).
Теперь задача построения передней стенки на втором этапе итераций сводится к отысканию гладкой кривой с максимумом угла наклона на оси луча (х = 0) и началом в точке (х = —Гр, г = 0). При этом глубина канала (г = гтах) определяется из условия, что на
этой кривой для температуры выполняется соотношение (5) при значении температуры на ней ТI _ Т8а4, где Т8а4 — температура кипения металла. В численной модели та-
I X— ¿1с
кая кривая находится приближенно. Здесь передняя стенка находится итерационно в виде кубического сплайна — путем нескольких решений краевой задачи. Итерации начинаются с некоторого нулевого приближения, которое затем поправляется в ходе итераций. В этом нулевом приближении задается положение передней стенки с глубиной, несколько большей, чем значение максимума по г изотермы кипения, полученной на первом этапе. В качестве начальной температуры на передней стенке берется температура кипения. Для поправок положения стенки имитируется описанный выше механизм саморегуляции ее положения. При этом полагается, что передняя стенка нижним концом касается окружности, форму которой имеет канал вблизи ее максимума по г. Радиус этой окружности г^ определяется из условия динамического равновесия свободной поверхности канала в окрестности его дна и выписан ниже. Точка касания находится либо в точке ж _ г^ — на границе луча лазера, либо (при малой мощности лазера) в точке, где угол наклона к оси ж строящегося сплайна становится близким к нулю, т. е. он не может быть продолжен дальше из требования того, что на нем Т _ Т8а4. В последнем случае максимум по г стенки канала находится в зоне, освещенной лучом. Тогда в этой точке сплайн заканчивается и касается упомянутой окружности.
Поскольку лазерная сварка существенно трехмерна, при наличии плазменного канала в процессе его движения имеет место перенос тепла конвекцией от передней стенки к задней. Частицы жидкого разогретого металла, находившиеся на передней стенке, в последующие моменты времени в результате переноса оказываются на поверхности задней стенки. В рамках тепловой плоской модели такой перенос тепла можно моделировать в записи граничного условия теплового баланса на задней стенке канала.
Для определения формы задней стенки условия (5) и сделанных допущений пока недостаточно. На втором этапе расчетов предлагается в качестве некоторого приближения к форме задней стенки принять отрезок прямой. Этот отрезок одним своим концом касается упомянутой окружности, а другой его конец совпадает с концом изотермы кипения на поверхности г _ 0.
Опустим подробности, связанные с тем, что уже на втором этапе положение концов стенок канала можно определить несколько точнее. Здесь только заметим, что вариация их положения влияет на расчет в небольшой их окрестности. Главным фактором, от которого зависит картина распределения температуры в расчетной области сварки, является тепловой процесс в окрестности луча лазера. Наклон передней стенки на оси луча и в ее окрестности определяет глубину канала.
В расчете на передней стенке в окрестности оси луча можно достаточно точно выполнить условие (5). Поскольку численное решение задачи без плазменного канала сходится на последовательности сеток с первым порядком, требуется, чтобы оно при наличии канала имело тот же порядок сходимости. Отсюда вытекает требование к точности выполнения условий теплового баланса на поверхности канала.
Дополнительные детали построения стенки проще изложить в связи с ее уточнением на третьем этапе решения задачи.
Условия на внутренних свободных границах служат для поправки их положения в ходе итераций, начиная со второго этапа. Условия выполняются с точностью порядка аппроксимации, когда сойдутся итерации. В данной работе удалось достаточно просто реализовать условия на внутренних границах. Реализованный способ состоит в следующем. В конце каждой итерации соотношение, аппроксимирующее условие в узловых точках внутрен-
них границ, разрешается относительно значения температуры в этой точке. В остальных точках берется значение температуры, определенное на данной итерации. В полученном дискретном поле температур описанным выше способом путем аппроксимации выстраивается уточненное приближенное положение этих границ. Затем определяются ближайшие к ним узлы сетки, как было описано выше. В узлах изотермы ликвидуса полагается T = 7]0, в узлах изотермы эвтектики — T = Te. Они служат краевыми условиями для температуры на координатных линиях x = const при численном решении задачи прогонками. В результате применения этого способа на следующей итерации внутри каждой подобласти решается задача с заданными значениями температуры на их границах. В пределе, когда сходятся итерации, на внутренних и внешних границах выполняются условия теплового баланса. В экспериментах на полученных численных решениях проверялась также точность выполнения условий на границах.
Уместно заметить, что совокупность первого и второго этапов представляет собой некоторую законченную приближенную численную модель лазерной сварки металлических пластин. Значение температуры, полученное в конце второго этапа расчетов, берется как начальное для третьего этапа.
5. Уточнение формы парового канала с использованием условия динамического равновесия
Для уточнения формы парового канала на третьем этапе кроме условия теплового баланса (5) предлагается использовать условие динамического равновесия на его поверхности. Оно записывается в виде
P (T ) = aKr + gpaz, (12)
т. е. на стенках канала сила давления паров уравновешивается силой поверхностного натяжения и гидростатического давления. В (12) a — коэффициент поверхностного натяжения жидкого металла (сплава); g — ускорение свободного падения; Kr — кривизна поверхности, определяемая соотношением
Kr = div | VZc = | . (13)
Vv1 + ivzc |2
Избыточное давление пара P является суммой статического давления Ps при поверхностном испарении, а также давления отдачи (реакции) Pr: P = Ps + Pr.
Давление насыщения смеси испарившихся компонентов сплава определим по следующей формуле [14]:
1 1 Р
Ps = £ 0.01 CjPjs = £ 0.01 Cj Pj eAj/T, (14)
j=0 j=0
где P0j, Aj, Bj — некоторые эмпирические константы. Для величины Pr примем линейную связь с давлением Ps: Pr = bPs.
Считая пар идеальным газом, для которого справедливо уравнение Клапейрона — Менделеева, получим выражение для плотности пара j-го компонента:
Pvj = RP2 /T, (15)
а для смеси —
^ (Т)_ ¿ 0.01 с ЯР ел/т, (16)
3—0
где — молярная масса; Я — газовая постоянная.
Соотношение (12) указывает на зависимость формы стенок парового канала от температуры. Здесь возникает сложная задача: необходимость решать уравнение теплопроводности (1), (2) с выписанными условиями теплового баланса на границах совместно с уравнением в частных производных второго порядка (12) на неизвестной границе парового канала. Предлагается упрощенный подход к вычислению поправки формы стенок канала.
Максимальному давлению на дне канала соответствует наибольшая кривизна. В силу достаточного поверхностного натяжения в этой точке поверхность канала близка к сфере с радиусом га, соизмеримым с радиусом луча.
Передняя стенка канала, построенная из требования выполнения условия (5) и касающаяся окружности радиуса г^, в данной модели больше не уточняется. Считается, что на ее форму в большей мере оказывает влияние испарение металла под действием лазера.
На третьем этапе уточняются приближенное распределение температуры и форма задней стенки канала. Отметим, что на втором этапе она взята в виде прямой, проходящей через конец полученной на первом этапе изотермы кипения (за лучом на границе г _ 0) и указанной окружности радиуса г^ на дне канала. Последовательно из (12) определяется кривизна задней стенки канала, а на следующей итерации в полученной расчетной области решается краевая задача с условием (5) на задней стенке. При этом учитывается только отраженный свет от лазера. Кроме того, кривизна задней стенки в плоскости у _ 0 мала по сравнению с ее кривизной в сечениях _ за исключением окрестности дна канала. По полученной кривизне подсчитывается радиус окружности в плоскости г7 _ ^2, касающийся задней стенки. Центр этой окружности полагается на оси, проходящей параллельно оси г через точку экстремума на дне канала. Во внутренних циклах, как и на втором этапе, вычисляется поле температуры в расчетной области последовательно во всех ее подобластях.
6. Оценки некоторых параметров парового канала
В расчетах в качестве материала пластин брали сплав АЛ2 (А1 +10% Б1).
Пусть г^ — радиус кривизны дна канала, который аппроксимируем частью сферы с этим радиусом. Тогда давление на дне
Р№)_ - + (17)
га
В этом случае, пренебрегая малой величиной гидростатического давления, получим
-
г« _ Рта • (18)
где Р(Та) _ Р3 + Рг _ (1 + Ъ)Р3 _ 1.55Т-1(0.9е33-294-37723Л4/т« + 0.1е43-584-63590-78/т), Ъ _ 0.55. Значение Та определяется из решения краевой задачи для температуры и составляет Та _ 2123 К _ 1850 °С. При Та _ 2123 К величина Р(Та) _ 3671.478 Н/м2. Следовательно,
0.57
га _-_ 0.15 мм.
а 3671.478
Таким образом, радиус кривизны дна канала составляет около 0.15мм, а давление на поверхности — около 0.037атм, т.е. около 4% от атмосферного давления. Для верхнего отверстия канала аналогично предыдущему имеем
1 + 1 = Рв
Гж Гх О '
Здесь гх — азимутальный, гх — аксиальный радиусы кривизны верхнего отверстия канала. Поскольку гх << г,г, характерный размер радиуса канала в этой области оценим как
О
гв
р (ТВ) •
Приняв Тв = 2003 К = 1730 °С, получим Рв = 1332.02 Н/м2, гв = 0.43 мм.
Если предположить, что на некотором расстоянии от оси луча температура на стенке канала достаточно равномерная, как и сила поверхностного натяжения и давления паров, то в некотором приближении окрестность дна канала можно считать близкой к сфере. Ее радиус определяется прежде всего величиной давления паров на дне канала и уравновешивающей его силой поверхностного натяжения. Полагаем приближенно, что задняя стенка в окрестности дна канала касается окружности с радиусом г а в плоскости у = 0, а ее сечение плоскостью г = 0 — окружности с радиусом г в и что центры этих окружностей находятся на оси, проходящей через вершину дна канала параллельно оси г. Таким образом определяется точка Мв с координатой г = 0 (у = 0) на задней стенке. В плоскости у = 0 заднюю стенку в некотором приближении можно аппроксимировать кривой второго порядка, касающейся в окрестности дна окружности радиуса г а и прямой г = 0 в точке Мв.
7. Численные эксперименты
Созданы программы решения рассматриваемой краевой задачи разностным методом, пригодные для моделирования процесса сварки в случае различных материалов и различных режимов. Проведены численные эксперименты на последовательности сеток с убывающими шагами прежде всего в расчетной области с фиксированной длиной. Наблюдались расчетные параметры задачи (распределение температуры, положение изотерм, границ парового канала и между фазами металла). Обнаружен первый порядок их сходимости при измельчении шагов сетки.
Описанный выше прием при выборе краевого условия на правой границе уменьшил его влияние на распределение температуры справа от ванны при длине расчетной области, при которой температура на правой границе не может быть близкой к температуре на бесконечности. Дополнительно для уменьшения влияния на численные результаты размера расчетной области в горизонтальном направлении эта величина варьировалась при одних и тех же значениях остальных параметров задачи и полученные результаты сравнивались. В данных расчетах удовлетворительным считалось, когда при увеличении размера области вдвое и независимом измельчении шагов сетки расчетные величины в областях парового канала и положение границ между фазами отличались не более чем на 3 %. При этом количество узлов расчетной сетки достигало полумиллиона. В этой задаче необходимый размер расчетной области, естественно, зависит прежде всего от скорости сварки. Чем она выше, тем больше несимметричность в картине изотерм относительно оси луча
Рис. 2. Поле температур в расчетной области: линия 1 — изотерма Т = 1207 К, 2 — линия ликидуса (Тю = 862 К), линия 3 — изотерма Т = 517 К.
лазера. Перед каналом (в направлении движения луча лазера) расстояние между изотермами меньше, чем между соответствующими им изотермами за каналом.
В данных расчетах физические параметры отвечают материалам, используемым в самолетостроении и других областях техники. Скорости сварки брались соответствующими реальному процессу сварки. Определялись температурные поля в изделии, положение внутренних границ между фазами материала изделия, форма и глубина парового канала. По результатам расчетов можно сделать прогноз о величинах зон, занимаемых в процессе сварки различными фазами, и установить, какая скорость движения луча лазера вдоль сварочного шва при заданной мощности лазера может обеспечить наличие достаточной по величине зоны жидкой фазы при отсутствии сквозного парового канала (пробоя изделия лучом лазера).
На рис. 2 приведены область и картина изотерм одного из расчетов (начало координат в левом верхнем углу расчетной области). За единицу длины принят радиус фокального пятна лазера. Для наглядности по осям х и г взяты разные масштабы. Длина расчетной области 660, координата х оси луча лазера равна 60. Здесь более низким температурам соответствуют более темные оттенки серого цвета. Находящаяся в левой части расчетной области маленькая темная подобласть в зоне высокой температуры, граница которой — изотерма 1 (через нее проходит луч лазера), соответствует паровому каналу.
8. О влиянии величины температуропроводности на размеры ванны с жидкой фазой
С практической точки зрения интересно решение проблемы о возможности лазерной сварки без образования парового канала. С одной стороны, его наличие приводит к определенным недостаткам этого процесса. Прежде всего, стенки канала как свободные поверхности жидкости слабо устойчивы. Колебания некоторых параметров процесса вызывают существенные осцилляции поверхности канала, что отрицательно влияет на однородность сварочного шва и его качество. Если сварку можно было бы осуществлять без возникновения парового канала, то процесс был бы более стабильным, более контролируемым, а шов — более однородным. Однако при небольшой мощности лазера ввиду недостатка энергии, подводимой в зону сварки, ванна с жидким металлом получается маленькой. Ее бывает достаточно для сварки только тонких пластин металла. При увеличении мощности используемого лазера из-за высокой концентрации энергии в окрестности оси луча
5,5
6,0
6,5
70
х
2-
4-
г
Рис. 3. Положение изотермы ликвидуса: при коэффициенте температуропроводности жидкой фазы, равном 100, 155 и 400 (кривые 1-3 соответственно).
металл начинает интенсивно испаряться. Гипотетически, если бы каким-либо способом можно было увеличить температуропроводность жидкой фазы, то тепло от окрестности оси луча лазера отводилось бы более интенсивно на периферию. Тогда верхний предел (критическое значение) мощности лазера, при котором в процессе сварки нет парового канала, можно было бы увеличить. При этом увеличились бы и размеры ванны. В реальных условиях с помощью внешнего физического воздействия (ультразвуком, электромагнитным полем и др.) можно вызвать микродвижения частиц жидкости в ванне и тем самым увеличить ее эффективную температуропроводность. При этом возникает вопрос: а какова зависимость критического значения мощности лазера от величины коэффициента температуропроводности жидкой фазы?
В данной работе исследовалось влияние коэффициента температуропроводности жидкой фазы на величину занятой ею ванны. При тех же значениях остальных параметров, при которых представлен результат расчета на рис. 2, на рис. 3 приведено положение изотермы ликвидуса (граница жидкой ванны в сечении у = 0) в той же расчетной области при трех различных значениях температуропроводности. Здесь по-прежнему ось луча лазера направлена вдоль оси г и имеет по оси х координату хс = 60 (единица длины — радиус фокального пятна). Получено, что глубина ванны в третьем приведенном на рисунке результате увеличилась примерно в 2.5 раза по сравнению с первым, т.е. толщина пластин, свариваемых в докритическом режиме, соответственно увеличилась в 2.5 раза. Однако следует отметить, что пластины все же достаточно тонкие.
Проведен также расчет при третьем значении температуропроводности жидкой фазы, когда описанные выше краевые условия на нижней поверхности пластин и верхней поверхности, за исключением свободной поверхности ванны, были заменены на условие равенства нулю потока тепла на них. При указанной скорости сварки эта замена привела к уменьшению критической мощности и увеличению линейных размеров ванны примерно на 10 %. Это говорит о том, что при данной скорости сварки потеря количества тепла через поверхности пластин в области ванны меньше, чем требуется для разогрева металла, поступающего туда вследствие движения со скоростью сварки.
Здесь уместно добавить, что в отдельных случаях возбуждение высокочастотных микродвижений частиц жидкости в потоке может привести к увеличению устойчивости некоторых течений жидкости. Известны эксперименты в аэродинамике, когда воздействие звука на поток увеличивает критическое число Рейнольдса, при котором ламинарное течение переходит в турбулентное. В данном случае не исключено, что при возбуждении в электропроводящей жидкости вихревых токов при определенных направлениях электродвижущей силы поверхность жидкой ванны будет более устойчива, чем при отсутствии токов. При этом дополнительная скорость движения частиц может увеличить эффективную теплопроводность в жидком металле.
Список литературы
[1] Лазерная и электронно-лучевая обработка материалов / Н.Н. Рыкалин, А.А. Углов, И.В. Зуев и др. М.: Машиностроение, 1985.
[2] Kroos J., Gratzke U., Simon G. Towards a self-constantent model of the keyhole in penetration laser beam welding //J. Phys. D: Appl. Phys. 1993. Vol. 26. P. 474-480.
[3] Судник В.А., Радаи Д., Дорофеев В.А. Компьютерное моделирование лазерно-лучевой сварки. Модель и верификация // Сварочное производство. 1997. № 1. С. 28-33.
[4] Semak V.V., Damkroger B., Kempka S. Temporal evolution of the temperature field in the beam interaction zone during material processing //J. Phys. D: Appl. Phys. 1999. Vol. 32. P. 1819-1825.
[5] Гладков Э.А., Гаврилов А.И., Малолетков А.В. и др. Динамическая нелинейная модель технологического процесса лазерной сварки с глубоким проплавлением // Сварочное производство. 2001. № 12. С. 17-24.
[6] Голубев В.С. Анализ моделей глубокого проплавления материалов лазерным излучением. Шатура, 1999. (Преп. / РАН. Ин-т проблем литья. № 83).
[7] Борисов В.Т. Теория двухфазной зоны металлического слитка. М.: Металлургия, 1987.
[8] Ораевский А.Н. Гауссовы пучки и оптические резонаторы // Тр. физ. ин-та им. П.Н. Лебедева. М.: Наука, 1988. Т. 187. C. 3-59.
[9] Кутателадзе С.С. Основы теплообмена. Новосибирск: Наука, 1970.
[10] Semak V., Matsunava A. The role of recoil pressure in energy balanse during laser materials processing // J. Phys. D: Appl. Phys. 1997. Vol. 23. P. 2541-2552.
[11] Басин А.С., Шишкин А.В. Получение кремниевых пластин для солнечной энергетики: Методы и технологии. Новосибирск: Ин-т теплофизики СО РАН, 2000.
[12] Зиновьев В.Е. Теплофизические свойства металлов при высоких температурах: Справочник. М.: Металлургия, 1989.
[13] Черепанов А.Н. Анализ подобия в процессах кристаллизации и структурообразования двойных сплавов // Металлы. 1988. № 3. C. 69-76.
[14] Черепанов А.Н., Шапеев В.П., Фомин В.М., Семин Л.Г. Численное моделирование теплофизических процессов при лазерно-лучевой сварке // Прикладная механика и техническая физика. 2006. № 5. C. 88-96.
Поступила в редакцию 13 марта 2006 г.
