Численное моделирование электронно-лучевой обработки материалов с учетом поверхностной активации и внутренних механических напряжений Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 15, № 3, 2010

Численное моделирование электронно-лучевой обработки материалов с учетом поверхностной активации и внутренних механических напряжений

А. Г. Князева

Институт физики прочности и материаловедения СО РАН, Томск, Россия

e-mail: anna@ispms.tsc.ru

A.B. Тян

Томский государственный университет, Россия e-mail: tav014@mail.tsu.ru

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

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

Введение

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

В последние десятилетия активно разрабатывались лучевые методы модификации поверхностных слоев материалов. Отметим, что электронные пучки обладают высокой плотностью энергии, достаточной для разогрева, оплавления и даже частичного распыления поверхностных слоев за весьма короткое время. Одно из характерных достоинств электронно-лучевой обработки — возможность применять импульсные пучки с контролируемыми временами длительностей импульсов и пауз. При этом выделяющаяся мощность может быть огромной (так, 1 Дж/см2 за 10 не соответствует удельной мощности 100 МВт/см2 [1]).

Известно, что на поверхности любого сплава в нормальных условиях формируется тонкий оксидный слой (например, 10-15 им для сплавов на основе никелида титана) [1, 2]. Модификация поверхностей сплавов электронным пучком приводит к насыщению прилежащих к поверхности слоев кислородом на глубину, почти в три раза превышающую толщину исходной оксидной пленки. Такая модификация позволяет осуществлять направленное перераспределение элементного состава приповерхностных областей, создавать более мелкую кристаллическую структуру и новые центры кристаллизации.

© ИВТ СО РАН, 2010.

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

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

Согласно экспериментальным данным [3, 4], в результате интенсивного внешнего воздействия на поверхность сплава в приповерхностной области образуется активированными слои ^ в котором наблюдается аномальное ускорение массопереноса. Эффект активации приповерхностных слоев объясняется тем^ что интенсивное внешнее воздеи-ствие в виде электронного пучка оказывает на атомы поверхности сплава не только мощное тепловое, но и механическое воздействие. В результате образующаяся волна механических возмущений приводит атомы в неравновесное состояние, которое характеризуется увеличением объема, приходящегося на один атом, возрастанием частоты колебаний, ослаблением межатомных и межмолекулярных связей, т. е. активацией поверхности. Фактически эффект активации поверхности является результатом "быстрых" процессов, протекающих за малые промежутки времени. В качестве дополнительного параметра, характеризующего степень отклонения системы от состояния равновесия, в работе [5] был введен параметр активации (параметр неравновесного состояния) П, определяемый по отклонению атомного объема от равновесного значения.

Как правило, в работах, посвященных моделированию процессов обработки материалов, не рассматривается роль механических напряжений, обусловленных процессами тепло- и массопереноса, в связи с чем может сложиться мнение об отсутствии напряженно-деформированного СОСТОЯНИЯ в диффузионной зоне в процессе обработки или, по крайней мере, о несущественном влиянии этого состояния на кинетику и механизм процесса обработки. Однако экспериментальные данные [3, 4, 6] указывают на то, что напряжения могут достигать величин, близких к пределу прочности материалов, и приводят к повреждению поверхности еще во время ее модификации. Первые экспериментальные данные о величине напряжений, образующихся при окислении металлов, были получены П.Д. Данковым и П.В. Чураевым [7]. Более поздние исследования показали, что эти напряжения достигают 1 ГПа [8, 9]. Подобные напряжения вызывают деформации, легко выявляемые методами металлографического и электронно-микроскопического анализов.

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

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

1. Математическая постановка задачи

Рассматривается одномерная модель процесса обработки образца сплава конечной длины Ь в декартовой системе координат. Математическая постановка задачи аналогично [10] включает в себя уравнения теплопроводности, диффузии кислорода и кинетики параметра активации с соответствующими граничными и начальными условиями.

Распространение тепла в сплаве с постоянными теплоемкостью с, плотностью р и коэффициентом теплопроводности Л описывается уравнением теплопроводности

дТ В 2Т

срт = ЛдХ2' (1)

справедливым в адсорбционном слое обрабатываемого сплава, где Т Ь — время, х — координата. Учитывая, что адсорбционный слой толщиной к является тонким по сравнению с характерным тепловым масштабом, т. е. выполняется условие к << лДГй, гд ,е кт = Л/(ср) — коэффициент температуропроводности, Ь* — характерное время, сопряженную двухслойную задачу теплопроводности можно свести к решению более простой задачи. Интегрирование уравнения теплопроводности по толщине адсорбционного слоя дает граничное условие на поверхности в виде

дТ дТ

х = 0 : -Л— = д(Ь) - £а(Т4 - Те4) - (ср)к-. (2)

Здесь

МЬ) = {

д(ь) = ЯоМЬУ,

1, (и + 1р){г - 1) < 1<и + (и + гр)(г - 1),

0, и + (и + Ьр)(г - 1) < Ь< (и + Ьр)г;

а — константа Стефана—Больцмана, Дж/(с-К4-см2); £ — степень черноты; Те — температура окружающей среды, К; ^ и Ьр — длительности импульса и паузы, с; г = 1,..., и, где и — число импульсов; до — плотность мощности теплового потока, которая представляет собой количество тепла ф, проходящее в единицу времени Ь и отнесенное к

-2

Выражение для потока массы для изотропной среды в случае приближения неидеального бинарного раствора имеет вид [11, 12]

3°х = -рО*ах911^сОх + аОх °х МТ °Х , (3)

где МОх — молярная масса кислорода; Я — универсальная газовая постоянная; аОх — линеиныи коэффициент концентрационного расширения для кислорода; В°х — коэффициент самодиффузии кислорода; дц — термодинамический множитель, в общем случае являющийся функцией состава. В первом приближении дц ~ 1. В этом случае

эффективный коэффициент диффузии кислорода О0х совпадает с коэффициентом самодиффузии: Б0х = О0х .

Компоненты первого инварианта тензора напряжений акк находятся из решения задачи о механическом равновесии.

В (3) учтены два механизма массопереноса: диффузия и перенос под действием градиента напряжений. В работах [5, 13] дополнительное ускорение массопереноса в зоне обработки связывается с уменьшением энергии активации диффузии Е0х вследствие эффекта активации поверхности. Математически эту зависимость можно представить в виде

Еох - IV

о

ох

О е*р(-Е°—[1) = До ехр(Уп),

где О0 = ехр(-Е0х/КТ) — коэффициент объемной диффузии в неактивированном материале; коэффициент 0 < = 1/&Т < 1 зависит от типа обрабатываемого вещества и легирующего элемента и определяется экспериментально, а величину 7 в первом приближении можно трактовать как энергию, необходимую для "активации" элемента, в расчете на один моль. Тогда для описания рассматриваемой диффузии одного компонента (кислорода) достаточно уравнения баланса массы

дС0х „ т /,ч

= ' ^ох, (4)

где 30х определяется выражением (3). Граничное условие на поверхности для уравнения диффузии кислорода (4) выводится так же, как для уравнения теплопроводности. Интегрируя уравнение диффузии по поверхностному слою к << \1О0хЬ* с учетом граничных условий (отсутствие потока массы на поверхности и идеальный контакт в точке х = к), получим граничное условие для уравнения (4) в виде

— П дС0х _ к дС0х а0хМ0хС0х д /кч

Х = 0: ^ = ~дГ + РЁТ дХакк- [Ь}

В простейшем приближении уравнение кинетики для параметра активации можно представить следующим образом [5]:

^ = кп фх(п)ф2(С0х), (6)

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

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

Ф1(П) = 1 - п-

Функция характеризует активацию поверхностного слоя по глубине:

{

г/ / \ I 1 Х/Ха •> х — ха •>

Ф " 1 (Х/Ха) Но, х > Ха,

ха

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

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

На торце образца при х = Ь примем условия отсутствия источников и стоков тепла и массы:

дТ

х = Ь : — = 0, Зох = 0. (7)

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

< = 0 : Сох ={£"Х>0: Т = Т" п = 0 (8)

Так как в модели рассматривается диффузия только одного элемента — кислорода, опустим в дальнейших рассуждениях индекс " Ох".

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

Чтобы найти воспользуемся подходом, описанным в [3, 6]. Обобщенные соотношения между компонентами тензоров напряжений и деформаций примем в виде

Оу = + (Хь£кк - Кт), (9)

где Аь, Ць — коэффициенты Ламе, К = Аь + /3 — изотермический модуль всестороннего сжатия,

т = 3 ат(Т - Т0) + ^ ак(Ск - Ско)

(к)

ат — линейный коэффициент теплового расширения, ак — линейные коэффициенты концентрационного расширения для компонента к, £кк = £и + £22 + £33-

Поскольку в данном варианте модели имеются диффундирующий кислород с концентрацией С = С и остальные элементы с концентрацией С2, запишем

С2 + С = 1,

следовательно,

т = 3 [ат(Т - То) + (а - а2)(С - Со)].

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

Окк = (3ХЬ + 2^ь)£кк - 3Кт = 3К(£кк - т), оц /А т

\ е т г

£И = 2^ - 5У I 6К^°кк - 3) : ИЛИ = + + :

где £еу — обратимые (упругие) деформации; £Т — термические деформации; £гу — деформации, связанные с необратимыми процессами переноса массы.

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

Полагая, что обрабатываемая пластина сплава не закреплена и свободна от действия внешних сил, а скорость переноса массы диффузией много меньше скорости распространения механических возмущений, придем к задаче о механическом равновесии, решение которой приводится в [3, 6] и в наших переменных имеет вид

£кк

1 - . „ , 1 + vw 2--(А1Х + А2) +

1 - V

1 - V 3'

акк

2Е I w

V А1Х + А2 - з

(10)

где V — коэффициент Пуассона, Е — модуль упругости, А1 и А2 постоянные пнте-грирования, являющиеся функциями времени:

ь ь

4 Г 2 Г

А1 = — - — wdx,

оо

ь ь

42 А2 = — wdx —— wxdx. 3^ Ь2 у

оо

В результате задача включает уравнения (1), (4), (6) с граничными условиями (2), (5), (7) и начальными условиями (8).

2. Безразмерные переменные и оценка параметров

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

т = г/г*, £ = х/х*, в = (Т - То)/(Т - То), Б^ = а^/а*,

где г* = Ьг, х* = Vкт г*, Т* = То + дох*/Х, а* = ЗКат(Т* - То) характерные масштабы. Концентрация С и степень активации п определению. Получаем постановку задачи в новых переменных:

дв = № дТ = д ё2 '

соответствующие безразмерные по

(П)

дС т д

дТ = Ье 61

(

О 1 +

2к1дшС\ дС 1 + вБг) д£

12Ье • к1ш д

(1 + в8г)Д2 д£

(АО С),

дП = К(1 - п)/(£/£а),

£ = 0:

дв В л дв

- -Ц = /г(т) [(1 + ввг)4 - 1] - 8—

' 2кгшдС\ дС 8 1 + вБг) ~д£

дС

дт

12к1шАС

Ье • В дт Д2(1 + ввг)

(12)

(13)

(14)

{ = А: f =0, * = 0,

-Г/ 2к\шдС\ dC 12hiuAC * = LeD[ ^ + TTMj dc " (1 + dSt)A2

(15)

г = 0: С

I

Со, С = 0, о С> 0,

в = 0, п = 0.

(16)

А

А

Обозначение А = дJ ю^С — J в котором ю = в(£,т)+ д(С ) — С 0)), введено

0 0

для удобства.

Безразмерный коэффициент диффузии

D-

ехр

eDn(St + 1) + St(e - 1)

в (i + est)

Выражение для объемных напряжений (10) в безразмерных переменных принимает вид

Skk = &kk/о* = 2hi

А А

(А - АО м+й - АО/

оо

wdC — w

(17)

Задача (11)—(16) содержит следующие безразмерные параметры: 8 = h/x* — толщина адсорбционного слоя; А = L/x * толщина образца; St = (T * — T0)/T0; B =

о , a — a2

Box * (T* — T0)o/A — параметр, характеризующий теплопотери; g = ——--- — от-

ат(T* - To)

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

атомных радиусов элементов. Комплекс h1 = --есть отношение скоростей рас-

Al +

пространения продольных и поперечных упругих волн. Так как продольные волны в твердых телах имеют скорость распространения больше, чем поперечные, то h1 < 1. Число Льюиса Le = D*/кт характеризует отношение скоростей тепловых и диффузионных процессов. Для твердых тел при T ^ T0 Le ^ 0 и при T > T0 Le < 1, так как скорость диффузии много меньше скорости распространения тепла. Комплекс aMo*

ш = ——--коэффициент связанности, кп — кп

— безразмерная скорость активации;

рЯТо

Са = ха/х* — глубина активации. Параметр = — безразмерная энергия, необ-

ходимая для активации поверхностного слоя. Параметр в = ЯТ*/Еэ — безразмерная температура активации диффузии.

3. Алгоритм численного решения задачи

В общем случае задача (11) (16) является связанной (т. е. в ней учитывается взаимодействие полей концентрации, температуры и напряжений) и для ее корректной численной

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

Задача теплопроводности решалась методом прогонки по безусловно устойчивой явно-неявной конечно-разностной схеме с ошибкой аппроксимации 0(йт1) + 0(й^2) [14]. Для реализации метода прогонки уравнение теплопроводности приводилось к виду

авг+1 — Ьвг + ¿вг-± + е = 0, (18)

где а = г; Ь =1 + 2г; с1 = г; е = вг; г = ¿т/а^2 — число Куранта. Температура вычислялась обычным способом [14]:

вг = «г+х^г+1 + вг+1, (19)

а ¿вг + е

где аг+1 = ---—, вг+1 = ,-1—, г = 1,п — 1- Граничные условия для прогоночных

Ь — ааг Ь — ааг

коэффициентов а1, /31 находились на основе разложения температуры в в точке г = 1

относительно точки г = 0 в ряд Тейлора со вторым порядком аппроксимации:

а1 = 2г^ • Z-1, 01 = + 25)во + 2!(т)ат — [(1 + ЗД)4 — 1] | •

где Z1 = 2гй^ + 25 + В точке г = п имеем

ХРп + V

вп

1 — апХ ''

= 2г = вп

ГДеХ = 2г + 1' ^ = 2г + 1-

Для решения задачи диффузии (12), (14)—(16) была использована двухшаговая двухслойная явно-неявная конечно-разностная схема, в которой расчет концентрации на текущем временном слое Сп+1 распадается на три этапа [15]. Сначала находится предварительное значение Сп+1, затемню той же схеме, но за два полушага вычисляется второе предварительное значение Сп+1. На заключительном этапе определяется окончательная величина с помощью "взвешивания", которое исключает главную часть погрешности:

Сп+1 _ 2 С п+1 _ С<п+1

Ст 2Ст Ст -

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

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

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

А + \А\ /С — С— М А — \А\/Сг+1 — С ^32С\ V ¿С +2 дС2) + 2 ^ ¿С 2 дС2)

2 \ ¿С А — \А\

+2 дС2 ) +

(

¿С

\А\^д 2С

\А\ А + \А\

Г< I 1___1 I

2¿С Сг+1 + 2^ Сг-1 + 2 дС2'

Для удобства преобразований в уравнении (12) введем обозначения:

№ ек1ш

( 2к1дшС\

1 + ОоБ^ ' (1 + ОоБг)^2

Старшую производную аппроксимируем следующим образом:

= V (г).

(20)

д /(С ) дС _ 1 дС1 (С,П) дС _ ¿С

/г+1 + / ■г Сг+1 — Сг /г + / г-1 Сг — Сг-1

2

¿С

2

¿С

В результате придем к разностной аппроксимации уравнения (12) в виде

Сг — С г ¿т

Ье

¿С

/ г+1 + / г Сг+1 — Сг / г-1 Сг — Сг-1

2

¿С

2

¿С

+

(

^1Ье/г \ Сг+1 — 2Сг + Сг-1

+ | —;--VDг\A\dС I -—-----„ „-Сг—

2

¿С2

(1 + е0Б1)в

— 2VD г

А — \А\С + \А\С А + А\С +1 + 1ССг —

(21)

Для уменьшения аппроксимационной вязкости при старшей производной вводим замену [15]

(

Ье/г

— VDг\A\dС

)

Ье/г

4VD г\А!С

2+

что корректно при условии

2VD г\А!С Ъе/г

о(1).

2

Это условие накладывает дополнительное требование на величину пространственного шага ¿С для диффузии. Тогда из (21) получаем

Сг — С г ат

Ге

¿С

/г+1 Сг+1 — Сг / г-1 Сг — Сг-1

2

¿С

2

¿С

+

Ге/г

Сг+1 — 2Сг + Сг—1

4УР г\Л\аС

2+

¿С2

2УЛ(БЬ +1)ев Р гПг (1 + воБ1)в

Сг - 2УРг

Л — \Л\С + \Л\С Л±Л!С

Далее приводим (22) к линеаризованной алгебраической форме (18), в которой

а= ¿=

Ь=

гЬе/г+1 2

гЬе/г-1 2

Ь гЬе/г • Z-1 — ЫСИгУ(Л — \Л\),

Ь гЬе/г • Z-1 + ЫСОгУ(Л + \Л\),

1 + а + й + Сг,

2УЛ(БЬ + 1)£пА пг ¿т в (1 + воБ1) '

(22)

Z2 = 2 +

24^шк^Л^С

Пг

дуг

Ж'

¿т

¿С2-

Д2(1 + воБг + 2шк1 Сг д)'

¿С г

кинетической задачах — разные, в то время как шаг по времени — одинаковый.

Граничные условия для прогоночных коэффициентов а1,в1 находятся так же, как в задаче теплопроводности:

1

а = гГе/о • Z',

о = С , 5 +1 е/о — 2ЛУРо)', в1 = ТС +2 Я/ ^

(23)

5 1 5^/—2ЛУВо) ¿т (Ге /о — 2ЛУ Р о)ЛУ Р о

^ = + о---+ 2ЛУDогdС--—-+

¿С 2 2 /о /о

ЛУРодт^п(БЬ + 1)по т , +--. „ ^--+ Ле/о.

в(1 + воБ1)

(24)

г=п ДЛЯ ЗЭДсХЧИ диффузии находятся так же,

как в задаче теплопроводности:

Сп

Х°п + V 1 — апХ;

Сп • Z- '

Z4 = 1 + 2гГе/п +

Х = 2Ле/п • Z- ' V

2ЛУРп£в (БЬ + 1^т 12шкЛ(2Г е/п^С + dт (Г е/п — 2ЛУР,п))

в (1 + воБЬ)

Д2(1 + воБЬ + 2шк1дСп)

е

_ / * _ /' _ /

/о _ дС} По _ дС} /п _ дС'

С помощью полиномов Лагранжа проведем одностороннюю аппроксимацию производных в (23)-(25) со вторым порядком точности [16]:

По _ ^ _ 2de +o(de

9fo 4fi - 3fo - f2 2\

f0 _ de _—щ—+o(dn

fL _ f _3fn - 4fn-1+fn-2+от.

de

2de

(26)

Формулы (23) справедливы во всей области определения параметров, кроме ситуа-2к1шд

ции, когда -- С _ —1 д < 0. В этом случае при старшей производной возникает

1 + воЬг

неопределенность вида (0 • то) и задача становится некорректной [17]. Для разрешения

2к1шд

проблемы прибегнем к следующему. Обозначим для удобства -- = ¡11 и вычислим

предел при g < 0:

ас

lim (1 + ßiC) —

1 + e0st

Иш

1 + ßC

ßi

с ^--

de Ci de/dC d2e/dC2

с=-i/ß!

( ö2C\

\ßi W)

C=-i/ß!

Таким образом, уравнение диффузии (12) при д < 0 в точке С _ —1/^1 преобразуется к виду

dC т д

dT _Le de

п д 2C

ßi ~W

д

- 2V— (ADC), ßi < 0,

совпадающему по форме с уравнением Кортевега—Де Вриза, которое описывает динамику уединенных волн (так называемых солитонов). Значение концентрации С _ —1/^1 является пределом растворимости диффундирующего элемента. Для данной одномерной однослойной задачи этот эффект может проявить себя только на поверхности, что приводит к необходимости иной аппроксимации граничных условий на поверхности, где будут справедливы уравнение

дС т Пе(Ы + 1)п' д2С т - д3С (Ьг + 1)п' дС

ж _ь** в(1+4 в? +Ье'" ° её— 2АК-ж+ьгС—А° вг

и условие

_ д2С _ 8 дС 6АС

С _0: ~д£2 _ - дТ

для которых аналогично предыдущему находим

ßi

ai _ 2AVgD0rde • Z-i,

^¡AVgdeS^n (St + 1)y'08 \ d3C

Co | -;--1—n/, , „ ^--1 -LeßiD0dr

\ ßiLe ß(1 + 9oSt)

de3

Z

i

(27)

бЛ2УдтдСД) _ ЛУд5дС

=-т^- + 2ЛУдД^С +-^ - 1-

Д2 ц\Ъе

£П(Бг + 1)^0(2ЛУд2Д2Ддт - №еЛ^ Ддт - 5дД2)

0 (1 + воБ1)Д2д '

(28)

д3С / ч

Производную — определим аналогично (26):

д 3С 1

= Щ?-^0 + 18С<1 - 24°2 + ЫСз - 3С4] + °(д?)'

Напомним, что в формулах численного счета диффузионной задачи (20)-(27) используется значение температуры на поверхности, поэтому уравнение (17) описывает термодиффузионные напряжения в микромасштабе диффузионной зоны. Если в решении задачи о механическом равновесии пренебречь концентрацией и учесть температуру по всему образцу, то после аналогичных операций обезразмеривания получим формулу для температурных напряжений в макромасштабе всего образца:

Бт = 2Нл

А А

(§ - ¿)/«+(Д - §) -,

9 = 9(т,С)' (29)

о о

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

1 а А

Б = Д У Бт + ДД У Бо

оо

где Бт определяются то формул е (29), а Б о — по формуле (17).

Уравнение кинетики параметра активации (13) можно проинтегрировать аналитически с помощью обычного метода разделения переменных. С учетом начального условия (16) получим

= кг, I (Шдт,

1 - п

-1п(1 - п) = кпI(С/Са) • т, п = 1 - ехр(-кп1 (С/Са) • т)' В силу особенностей рассматриваемой задачи, в частности сильной нелинейности уравнений, для выбора расчетных шагов по координате и времени необходимо проверить численное решение на сходимость, а также провести сравнение численного решения с точными аналитическими решениями простейших задач диффузии и теплопроводности. Поскольку в данной модели отсутствуют внутренние источники и стоки массы, при определении сходимости численного решения учитывали также выполнение интегрального закона сохранения массы

А

С(т, 0) • 5 + У С(т, С)дС = Со • 5, (30)

о

где Со — начальное значение концентрации в поверхностном адсорбционном слое 5.

Для выбора оптимальных пространственно-временных шагов на сетках диффузион-но-кпнетпческой и тепловой задач построены зависимости значений температуры и концентрации на заданной глубине от квадрата пространственного шага при различных шагах по времени для одного импульса к заданному моменту времени наблюдения. Из анализа построенных графиков выбраны шаги по пространству для температуры йе = 0.04, для концентр ации йе = 4 • 10-5, по времен и — йт = 3 • 10-2. При этом для температуры и концентрации относительные погрешности не превышали 1 %. Относительная погрешность выполнения закона сохранения массы (30) не превышала 5 %.

Если пренебречь в граничном условии (14) потерями тепла и принять длительность импульса бесконечной (источник непрерывного облучения), то для уравнения теплопроводности (11) с начальными и граничными условиями (14)—(16) можно получить точное аналитическое решение методом изображений Лапласа [18]:

0(£, т) = 2*- (5 + еМсШ + ^(Ш + ^Г) •

V п Vт \2л/т 5 )

Если в уравнении диффузии (12) с начальными и граничными условиями (14)—(16) опустить эффекты активации и напряжений, то для данной задачи при условии постоянства коэффициента диффузии также можно получить точное аналитическое решение:

г к ) г (+ Ье-Р т) , (Шё£т + е ^

г(е,т) = с»ехр I—*—) • ег,с (—т~ + шж) •

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

4. Анализ результатов

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

На рис. 1-4 представлены примеры численного исследования электронно-лучевой обработки на основе предложенной модели. Различия в распределениях концентраций кислорода и напряжений в диффузионной зоне для связанной и несвязанной задач иллюстрирует рис. 1. Активация материала приводит к существенному увеличению глубины и скорости диффузии. В частности, при данном наборе параметров в случае активированной диффузии величина остаточной концентрации на поверхности к моменту времени т = 0.8 почти на 50% меньше, чем в случае неактивированной диффузии. Из рисунка также видно, что связанный характер процессов переноса и деформирования приводит к существенному изменению концентрации кислорода в диффузионной зоне как для активированной диффузии (штриховые линии и символы о), так и для неактивированной (сплошные линии и символы А). В квазистатической постановке напряжения полностью определяются концентрацией и температурой, поэтому характер распределения напряжений в диффузионной зоне при варьировании ео и ш аналогичен характеру распределения концентрации (см. рис. 1).

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

Ст

ности, на котором температура достигает некоторого заданного значения 9*. Зависимости Со и СТ от времени, определенные таким образом, показаны на рис. 2. Из рис. 2, а следует, что глубина диффузии растет до тех пор, пока температура при данном наборе параметров в течение действия импульса достаточна для диффузии, т. е. до момента т = т = 1. Характер поведения Ст зависит от 9* (рис. 2, б): при 9* = 0-1 наблюдается постепенное увеличение глубины прогрева, при 9* = 0-7 вначале наблюдается увеличение Ст

ся практически с начала процесса обработки, то во втором она "фиксируется" гораздо позже.

Рис. 1. Распределения концентраций кислорода (а) и напряжений (б) в диффузионной зоне в моменты времени т = 0.8 (1), 0.9 (2), 1 (3). Остальные параметры фиксированы: Ье = 0.01, 5 = 0.005, А = 15, £а = 8, Бг = 3, В = 10"5, Л-1 = 0.2, д = -10, кп = 0.1, в = 0.06 для одного импульса

Рис. 2. Зависимость глубины диффузии (а) при е = 0.1 (1), 0.3 (2), 0.9 (3) и глубины прогрева (б) при В = 10_3 (!), 10"2 (2), 10"1 (3), 1 (4) от времени; остальные параметры те же, что к рис. 1

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

На рис. 4 приведены результаты электронно-лучевой обработки поверхности двумя последовательными импульсами при тг = тр = 1. Видно, что поведение термических напряжений непосредственно зависит от режима облучения. Так, повышение температуры сопровождается ростом термических напряжений (рис. 4, а). Толщина адсорби-5

ясняется тонкостью этого слоя и его высокой теплопроводностью. Из рис. 4, б следует

02 04 М 08 L0

Рис. 3. Зависимость глубины диффузии от коэффициента связанности ш к моменту времени т = 5; ö = 0.006, Le = 0.001 (1), 0.005 (2), 0.007 (3), 0.01 (4); остальные параметры те же, что к рис. 2, 3

а б

Рис. 4. Изменение во времени: а — температура в (сплошные линии) и термические напряжения Бт (пунктир), б— средние напряжения >5; 5 = 0.001 (1), 0.01 (2), 0.1 (3) и обработке поверхности двумя последовательными импульсами; Ье = 0.01, А = 20, = 12, БЬ = 3, В = 10-3, Н\ = 0.2, д = -5, кп = 0.05, в = 0.05, ев = 0.1, и = 0.5

немонотонный характер кривых средних по образцу напряжений, что в окрестности т = 3 связано с различным поведением полей концентрационных и термических напряжений во времени.

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

Список литературы

[1] мейснер Л.Л. Механические и физико-химические свойства сплавов на основе никелида титана с тонкими поверхностными слоями, модифицированными потоками заряженных частиц: Дис. ... д-ра физ.-мат. наук. Томск, ИФПМ СО РАН, 2004. 546 с.

[2] Лотков А.И., мейснер Л.Л., Гришков В.Н. Сплавы на основе никелида титана: Ион-но-лучевая, плазменная и химическая модификация поверхности // Физика металлов и металловедение. 2005. Т. 99, № 5. С. 66-78.

[3] Еремеев B.C. Диффузия и напряжения. М.: Энергоатомиздат, 1984. 180 с.

[4] Гегузин Я.Е. Диффузионная зона. М.: Наука, 1979. 344 с.

[5] Князева А.Г., Псахье С.Г. Диффузия элементов в активированном поверхностном слое // Физическая мезомеханика. 2006. Т. 9, № 2. С. 49-54.

[6] Боли Б., УэЙНЕР Дж. Теория температурных напряжений. М.: Мир, 1964. 517 с.

[7] Данков П.Д., Чураев П.В. Эффект деформации поверхностного слоя металла при окислении // Докл. АН СССР. 1950. Т. 73, № 6. С. 1221-1224.

[8] Кубашевский О., Гопкинс Б. Окисление металлов и сплавов: Пер. с англ. В.А. Алексеева. М.: Металлургия, 1965. 428 с.

[9] Bradhurst D.H., Heuer P.M. The influence of oxide stress on the breakaway oxidation of ZircaIloy-2 //J. Nucl. Materials. 1970. Vol. 37, No. 1. P. 35-47.

[lOj Тян A.B., Князева А.Г., Псахье С.Г. Нелинейные эффекты в поверхностном слое никелида титана в условиях его неравновесной активации импульсным электронным пучком // Изв. вузов. Физика. 2007. Т. 50, № 3. С. 8-16.

[11] Князева А.Г. Диффузия и реология в локально-равновесной термодинамике // Математическое моделирование систем и процессов: Сб. науч. тр. / Под ред. П.В. Трусова. Пермь: Изд-во Пермского РТУ, 2005. № 13. С. 45-60.

[12] Князева А.Г. Перекрестные эффекты в твердых средах с диффузией // ПМТФ. 2003. Т. 44, № 3. С. 85-99.

[13] Князева А.Г., Псахье С.Г. Моделирование неравновесной диффузии, сопровождаемой внутренними напряжениями // Физическая мезомеханика. Спец. выпуск. 2005. Т. 9, № 8. С. 41-44.

[14] Самарский A.A. Введение в численные методы. М.: Наука, 1982. 271 с.

[15] Пасконов В.М., Полежаев В.И., Чудов Л.А. Численное моделирование процессов тепло- и массопереноса. М.: Наука, 1984. 288 с.

[16] миньков С.Л., миньков Л.Л. Основы численных методов: Уч. пособие. Томск: Изд-во НТЛ, 2006. 260 с.

[17] Князева А.Г. Режимы развития из начального зародыша твердофазной реакции, лимитируемой диффузией // ФГВ. 1996. Т. 32, № 4. С. 72-76.

[18] Лыков A.B. Теория теплопроводности. М.: Высшая шк., 1967. 600 с.

Поступила в редакцию 18 мая 2009 г., с доработки — 17 декабря 2009 г.