CC BY
49
2012 Механика № 3
УДК 531+669.539+536
Е.С. Ильина1, В.Н. Демидов1,2, А.Г. Князева1,2
1 Национальный исследовательский Томский политехнический университет, Томск, Россия 2Институт физики прочности и материаловедения СО РАН, Томск, Россия
ОСОБЕННОСТИ МОДЕЛИРОВАНИЯ ДИФФУЗИОННЫХ ПРОЦЕССОВ В УПРУГОМ ТЕЛЕ ПРИ ЕГО ПОВЕРХНОСТНОЙ МОДИФИКАЦИИ ЧАСТИЦАМИ
В работе предложена изотермическая динамическая модель начальной стадии процесса перемешивания частиц в поверхностном слое материала в условиях поверхностной обработки. Модель учитывает диффузию ионов, конечность времени релаксации потока массы, появление напряжений вследствие изменения состава поверхностного слоя и явление переноса компонента под действием градиента напряжений. Описан выбор безразмерных переменных и физический смысл получаемых при этом безразмерных комплексов. На основе анализа известных частных уравнений и классических разностных методов их численного решения выбран метод численной реализации связанной модели. Продемонстрированы свойства разностных схем и качественные особенности решений диффузионного, телеграфного и волнового уравнений. Выявлены качественные особенности концентрационной и упругой волн, взаимодействующих между собой. Показано, что смена качественного характера распределения деформаций в упругой волне связана с конечностью скорости распространения диффузионной волны.
Ключевые слова: перемешивание, диффузия, концентрационная волна, упругая волна, связанная модель, поверхностная обработка, разностные схемы.
E.S. Il’ina1, V.N. Demidov1,2, A.G. Knyazeva1,2
1National Research Tomsk Polytechnic University, Tomsk, Russian Federation institute of Strength Physics and Material Science, Tomsk, Russian Federation
THE MODELING FEATURES OF DIFFUSION PROCESSES IN ELASTIC BODY UNDER PARTICLES SURFACE TREATMENT
Isothermal dynamic model is suggested for the initial stage of the intermixing process in surface layer of material under surface treatment condition. He model takes into account the particle diffusion, the finite of mass flux relaxation time; the stress appearance due to composition change of surface layer and mass transfer phenomenon under stress gradient action. The choice of dimensionless variable and physical sense of obtained dimensionless complexes are described. The method of numerical realization of coupling model is selected on the base of the analysis of known particular equations and classical ways of their numerical solution. The features of difference schemes and qualitative distinctions of the diffusion, telegraph and wave equation solutions are demonstrated. The qualitative regularities are revealed in the concentration and elastic waves interrelating between each other. It is shown that the
change of the qualitative strain distribution in elastic wave connects with the finiteness of the propagation velocity of the diffusion wave.
Keywords: intermixing, diffusion, concentration wave, elastic wave, coupling model, surface treatment, difference scheme.
Введение
Взаимодействие мощных импульсных пучков заряженных частиц с твердым телом активно исследуется на протяжении последних трех десятилетий. Первые публикации на эту тему появились практически сразу после создания импульсных ускорителей и привлекли внимание специалистов. В них было показано, что при переходе от непрерывного к высокоинтенсивному импульсному режиму облучения происходит качественное изменение природы радиационно-стимулированных процессов. Данное обстоятельство помимо чисто научных задач породило значительные надежды на возможность технологического применения импульсных ускорителей заряженных частиц для модификации физико-химических свойств материалов и изделий [1]. В настоящее время различные научные школы активно работают в этой области. Различные виды поверхностной обработки различаются как длительностью и частотой импульсов, энергетическими характеристиками частиц, составом модифицирующих потоков и плазмы, так и техническими особенностями.
Поверхностная модификация сопровождается различными физическими и химическими факторами, которые непосредственно оказывают влияние друг на друга. Например, внедряемый в поверхностный слой поток ионов способствует нагреву, изменению состава вследствие перемешивания и образования соединений и фаз, формированию поля напряжений и т.п. Экспериментально исследовать роль каждого фактора в отдельности невозможно. Поэтому большую роль в их изучении в данном случае играет компьютерное моделирование.
При построении математических моделей разные авторы уделяют внимание различным сторонам процесса поверхностной обработки, как правило, не придавая значения совместному протеканию разных явлений. Например, в [2-4] моделируются неравновесные процессы плавления и кристаллизации с образованием Кнудсеновского слоя в приповерхностной зоне; обсуждаются вопросы корректной формулировки граничных условий во фронте испарения; исследуется влияние формы импульса на динамику процессов плавления и испарения. В ряде работ
этих же авторов и в [5-6] изучаются упругие и упругопластические волны, возникающие вследствие быстрого ввода энергии в вещество в условиях облучения. Но процесс внедрения ионов в этих работах остается без внимания. В отдельных публикациях рассчитываются пробеги частиц разной массы и энергии в средах различной плотности и структуры, но без учета сопутствующих явлений эти работы не описывают реальный эксперимент.
Цель данной работы - исследование начальной стадии диффузионных процессов в упругом теле при его поверхностной обработке частицами в изотермическом приближении.
1. Общие соотношения
В первом приближении процесс внедрения ионов можно считать изотермическим, если рассматривается низкоэнергетический поток (в этом случае изменения температуры невелики [1, 7]), а возникающие в зоне воздействия механические напряжения - упругими. Тогда для описания взаимодействия концентрационных и механических полей в условиях обработки однокомпонентного материала потоком ионов, состоящим из частиц другого сорта, требуются уравнение неразрывности, уравнения движения и уравнение баланса массы для компонента (внедряемых ионов):
^ + РУ-V = 0,
а?
¿с ^ т
р— = -V- з, а? ¿V V Р а
где р - плотность среды; С - массовая концентрация внедряемого компонента; I - поток массы этого компонента; а - тензор напряжений; V - вектор среднемассовой скорости. Для замыкания системы уравнений требуется выписать определяющие соотношения. В данном случае это относится к потоку массы и соотношениям между компонентами тензоров напряжений и деформаций. Если деформации малы, а напряжения можно считать упругими, то для изотропной среды имеем определяющие соотношения теории массоупругости [8]
ау = 2И£у +5У 1Л£кк - Км\, (1)
где екк = еп + е22 + е33; = 3Да(С - С0); Х,ц - коэффициенты Ламе;
2
К = X + з ц - упругий модуль всестороннего сжатия; Да - разность
коэффициентов концентрационного расширения внедряемого элемента а и основного материала а0. Компоненты тензора малых деформаций определяются соотношениями Коши:
dxj dxl
В приближении неидеального раствора для диффузии по механизму внедрения выражение для потока массы имеет вид [8]
J = -Df (С)VC + BCVatt -1, dji-, (2)
D0 mj \ . .
где B =----------(a-a0) - коэффициент переноса под действием напря-
RTp
жений; aкк =Оц +а22 + а33; Do - коэффициент самодиффузии, R -универсальная газовая постоянная, m - молярная масса; T - температура; tr - время релаксации потока массы; вид функции f (С) зависит от структуры образующегося в процессе имплантации «раствора».
2. Упрощенная одномерная модель Так как деформации малы, то
V’ V = dîtSkk и P = p0exp(-8») * Po •
Кроме того, предположим, что малы не только деформации, но и скорости, и ускорения. Тогда можем записать
dv
P d7 *po
d V — + vVv
Vdt J
d 2u
Po
dt
2 •
При условии, что поток равномерно распределен по обрабатываемой поверхности, придем к одномерной задаче.
В приближении одноосного нагружения, а11 Ф 0, а22 = а33 = 0, из (1) находим
а кк - а11 - Е
.811 - т
V 3.
Тогда
х -- д г (с + вс Х.
ОХ ОХ ОХ
(3)
(4)
В приближении одноосной деформации, еп Ф 0, в22 -в33 - 0, из (1) следует
г
а11 - 3К
X + 2ц w
------- £11----------
V 3К 11 3 У
) 3К 3К 4
акк - 3К (£11 - ^-^^а11 -^^-Цw.
X + 2ц X + 2ц 3
Подставляя последнее в выражение для потока массы (2), найдем
А>. Г (с)-
12 Кцв
(а-а о )с
X + 2ц
Если ввести обозначения
X + 2ц
Ос
Ох
+ в
3К с Оа
11
£11--
еп, в'-в-
X + 2ц Ох
3К
О1
ОХ '
/ '(с)- / (с)
3К 11 X + 2ц
12 КцВ (а - а 0)
+ -
X + 2ц
то оба варианта (одноосное напряженное и одноосное деформированное состояния) с математической точки зрения окажутся эквивалентными, что удобно для подробного параметрического исследования связанной задачи.
В результате одномерная модель примет вид
ос о
Ох
Ох
Х --До/(сус + в>с£ - Хг Х
Ох Ох ОХ
О 2и _ Оа ^ОХ2 Ох ’
(5)
(6) (7)
а- Е '(е-Да(с - с0)),
(8)
где ось (Ох) соответствует направлению облучения; а, в, и - компоненты тензора напряжений, деформаций и вектора перемещений в том же направлении. В (6) под коэффициентом в0 подразумевается в или В'; в (8) под величиной Е' - соответственно Е или К, а под величиной в - либо ВЦ, либо в11.
Дифференцируя (7) по координате, а (8) - дважды по времени, после несложных преобразований запишем уравнение движения в напряжениях:
р О 2 а
+ рДа
О2с О2а
Ох
(9)
Дифференцируя (5) по времени, (6) - по координате и избавляясь от смешанных производных, получим
О2с Ос О
Хг —-—I--------- —
Ох 2 Ох Ох
А / (с )Ос
Ох
В —
Ох
сОа
Ох
(10)
Граничные условия к (9), (10) запишем в виде х - 0: I - т0ф1 (х); а-а0ф2 (х),
с - 0, в-0,
(11)
(12)
где а0 «т^р; функции ф1 и ф2 определяются характером внешнего воздействия (облучения) и взаимодействием частиц с поверхностью образца. В простейшем приближении
Ф1 (х )-Ф2 (х ) •
Вследствие различия характерных физических масштабов разных явлений действие потока массы и механических напряжений могут быть несогласованны.
В начальный момент времени имеем
t = 0: C = 0, ст = 0, — = 0, — = 0.
д1 дt
(13)
Мы получили простейшую связанную нелинейную задачу массо-упругости.
3. Безразмерные переменные
В общем случае выбор масштабов - задача нетривиальная, требующая внимательного подхода. Остановимся на этой процедуре подробнее. Введем обозначения для новых (безразмерных) переменных:
„а т „ х в
S =—; т = —; 5 = —; e = —, с * ,* х* в *
(14)
где величины с индексом «*» пока не определены. Подставляя (14) в уравнение (10), получаем
8 20 'г ,2 8т 2
8С
,*8т
А) 8
2 х 85
/ (с )§
85 .
Бс* 8
х*2 85
с 85 8^
или
Примем
г Л 2
А0,* 8
2 х* 85
/ (с )8С
1 ;8^
4Бс, 8
х* 85
с 85 . в5.
Do ,* х2
= 1,
(15)
тогда масштабы для координаты и времени будут связаны соотношением
х* =
Подставляя (14) в определяющее соотношение для напряжений (8), получаем
Е 'в*
с*
е -
(а-а 0 )
(с - С0)
Примем за масштаб для деформации в* = а0. Тогда в качестве масштаба для напряжений удобно взять
2
*
а* = Е'а 0 .
В результате найдем
где
5 - е -у(С - Со),
(а-а о)
у = .
а
о
Множитель при втором слагаемом в уравнении (15) теперь принимает вид
,*Бс* Бс* Б
= — Е а.
х*2
А А
Используя определение Б, получим
Б А0та0Е'а0 ^ \ тЕ'а2
—Е 'а 0 = Ал
ЯГра 0 А0
(а-а 0) = -
ДГр
у = юу .
где ю =
тЕ 'а 2 ДГр
В результате уравнение (15) запишем в виде
д2С дС д
г а_2
Эх Эх д%
г (С)—
юу
_д
д5
с д5 . д£
Уравнение движения (9) преобразуется следующим образом:
р а* д 5
~еЦ эХ2
+ Р-
а
,(а-а0)д2С = а* д25
а 0^*
дх 2 х2 д£, 2
х2
Домножим уравнение на — .
а*
Тогда
р х2 д25 а0(а-а0)х2 д2С = д25
Е' г*2 дх 2
2
а 0^* а*
дх
2
д5
2
Подставляя сюда полученные выше выражения для масштабов, получим
Р А0
Е' и
д_Б 8т2
- + у-
8 2сЛ
8 2 Б 85 2
Поскольку мы можем произвольно выбрать масштаб времени, примем за единицу множитель, стоящий перед скобкой. Тогда
,* = — А0. Это даст нам второе уравнение Е
8 2 Б
8т2
- + у-
8 2с 8 2Б
8т2 852
Выбранный масштаб времени ,* - это время, за которое упругая волна пройдет расстояние, равное эффективной ширине диффузионной зоны.
Аналогичным образом обезразмеривались граничные условия.
В итоге после перехода к переменным (14) задача (9)-( 13) примет вид
82с 8с
2 ■ + ■ 8т2 8т
8
85
/ (с )8с
85
-юу
8
85
с
85
82 Б 8 2с = 82 Б
8т2 +У8т2 = 85^
Б = е-у(с - с0),
5 = 0: 3 = -/(с)—+ сюу— -тг —, 3 = рф(т), ^85 85 г 8т п ’
Б = Б0ф(т),
5^да: с = 0, Б = 0,
8с 8Б
т = 0: с = 0, Б = 0, — = 0, — = 0,
8т 8т
где ю, у, тг, р - параметры модели:
(16)
(17)
(18)
ю =
3тЕа0 ЯТ р
а
3,
Ап
г
*
Заметим, что в предложенной упрощенной модели нет ограничений на величину концентраций. Поэтому следует ожидать, что не при всех возможных значениях параметров модели решение имеет физический смысл.
В качестве исследуемых материалов были выбраны алюминий (имплантируемый материал) и титан (материал основы). Их свойства представлены в таблице [9, 10].
В итоге были определены интервалы значений ю, у, тг, Р :
тг е[10 -520], юе [0,0 -104 ],
|у|е[0,02 - 0,6],
Р е [0,004-10-1 -0,003-10-1 ].
Свойства исследуемых материалов
Свойство Материал
А1 Т1
Плотность, р, кг/м3 2,69-103 4,54 -103
Молярный объем, Ут, м3/моль 8,99-10-6 -10,0-10-6 10,6-10-6
Молярная масса, т, кг/моль 26,98 -10-3 47,88 -10-3
Модуль упругости, Е, ГПа 70-79 120-135
Коэффициент концентрационного расширения, у 0,458-0,485 0,556-0,514
Таким образом, задача имеет особенность, отличающую ее от аналогичных задач термоупругости. Параметр тг может меняться в
широких пределах - от малых до достаточно больших величин. При исследовании медленных процессов диффузии в равновесных условиях и в однородных или почти однородных средах роль времени релаксации незначительна. При изучении процесса переноса массы в интервале времени, характерном для распространения механических волн, конечность времени релаксации становится важной.
В принципе в этом случае в качестве масштаба времени можно выбрать и время релаксации t* = (г. Тогда вместо (15) получим
д2С дС Б0 д
—г + — -1,-0—
дх дх
_0
г х*2
х* д%
/ (с )дС
Во д
г х2 д^
с дБ
Пространственным масштабом в этом случае будет
х* — д//,^0 ,
расстояние, которое пройдет диффузионная волна за время . Уравнение движения примет вид
д 2 Б
дх2
+ у
д2С Е'гг д2Б
дх2
РА) д^2
То есть в этом случае уравнение диффузии будет без особенностей, а в уравнении движения появится большой параметр, который создаст дополнительные проблемы. Множитель перед производной в правой части есть квадрат отношения скорости упругой волны и скорости диффузионной волны:
Е 'гг РБ0
\2
4. Выбор разностной схемы для частных вариантов модели
Для численной реализации модели удобно перейти к деформациям, оставляя в каждом из уравнений производные по времени лишь для одной искомой величины. Используя (18), вместо (16), (17) получим
д2С дС
—V +--------
дх дх
д 2е дх2 '
д д5
д2е - дС д^2 у д^2 ’
* (с )!
юу
д5
С ^
д^.
Выражение для потока массы в этом случае примет вид
( ^дС де д1
—- * (с VСюуд^-х' * •
(19)
(20)
(21)
В (20), (21) * (С) — /(С) + юу2С - заданная функция концентраций.
V
е
,
За кажущейся «простотой» модели скрываются вычислительные трудности. В основном это связано с тем, что характерные пространственные и временные масштабы диффузионных и механических процессов - различны. При независимом протекании такие процессы описываются соответственно параболическим и гиперболическим уравнениями, для решения которых корректно применять специальные методы.
Устремляя поочередно параметры к нулю, мы приходим к известным уравнениям, методы решения которых хорошо разработаны. Это уравнение диффузии
телеграфное уравнение
и волновое уравнение
Вспомним качественные особенности решений простейших задач для этих уравнений.
Решения чисто диффузионной задачи при заданном постоянном потоке (граничные условия второго рода)
дС
-д^ — Р, 5 — 0 (25)
д5
и заданной постоянной концентрации на левой границе (граничное условие первого рода)
С — Сб , (26)
представлены на рис. 1.
дС _ д2С дх д5
2
д 2С дС — д 2С г дх2 дх д^2
д2е — д2е дх2 д^2
(22)
(23)
(24)
Точные аналитические решения этих простейших задач известны [11].
Численное решение осуществляли по неявной и явной разностным схемам. В любом случае численное решение (сплошные линии) и аналитическое решение (точки) практически совпадают.
Рис. 1. Сравнение численного (квадратики) и аналитического (сплошная линия) решений диффузионных уравнений с граничными условиями первого (а) и второго (б) рода; моменты времени Т : 1 - 0,1; 2 - 0,5; 3 - 1,0
Для численного решения телеграфного уравнения (23) использовали такие же типы разностных схем, что и в случае уравнения диффузии. Неявная разностная схема
Сп+і 2Єп + Сп~і Сп+і Сп
г Ат2 Ат
. /'тп+і /'тп+і . /'тп+і /'тп+і
Зі+і + Зі сі+і - сі Зі + 8і-і сі - 4-і
2
А5
2
А5
(27)
где і - номер пространственной точки; п - номер временного слоя. Явная схема «крест»
Сп+і - 2Сп + Сп-і Сп+і - С
п-і
іп+і
Ат2
Ат
л_
А5
~ І ~ /-<п Ґ~іп _ І _ Ґ~іп Ґ~іп
§і+і + Зі Сі+і - Сі Зі + Зі-і Сі - Сі-і
2
А5
2
А5
(28)
где I - номер пространственной точки; п - номер временного слоя.
К каждому заданному значению времени вещество проникает на конечное расстояние от поверхности, что принципиально отличает решение телеграфного уравнения от чисто диффузионных, где значение
г
концентрации стремится к начальному значению лишь асимптотически. В этом случае мы можем говорить о глубине проникновения примеси.
В задаче с граничным условием первого рода явная схема приводит к появлению осцилляций в местах разрыва решения, что показано на рис. 2 для трех моментов времени.
0.0 0,5 1,0 1,5 £ 10-2 0,0 1,0 2,0 3,0 ^ 1()-2 0,0 1,0 2,0 3,0 4,0 ^ ц)-?
а б в
Рис. 2. Пример решения первой краевой задачи для телеграфного уравнения,
Св = 0,2: т = 0,43-Ю"4 (а); т = 0,26-і0"3 (б); т = 0,43-і0"3 (в)
Осцилляции не имеют физического смысла, а есть особенность численной реализации с использованием явной схемы, которая условно устойчива, если соблюдено условие Куранта-Фридриха-Леви [і2]. Явная схема хорошо передает скорость распространения диффузионной волны и глубину проникновения примеси.
Неявная разностная схема имеет другую особенность (рис. 3): в окрестности разрыва концентрации наблюдается «диффузионное размытие» численного решения. И это «численная диффузия», а не диффузия физического происхождения, которая проявляется уже для больших значений времени. С уменьшением шага по пространству и соответственно по времени описанный эффект размытия заметно уменьшается.
Решения волнового уравнения исследованы для различных форм импульса. При тестировании задачи были использованы три типа разностных схем [і2]: неявная, неявная с весами и явная схема «крест».
Схема «крест»
(п+і О п . п-А і / п 0^1 п \
^ Є - 2е, + е ^^Не+і - + е-\
0,5 1,0 1,5 £
б
Рис. 3. Сравнение численного (квадратики) и аналитического (сплошная линия) решений телеграфного уравнения с граничными условиями первого рода (26) (а) и решение телеграфного уравнения с граничными условиями второго рода (25) (б); моменты времени г: 1 - 0,06; 2 - 0,3; 3 - 0,6, 4 - 0,9; 5 - і,2; 6 - і,5, 7- 6,0; 8 - 7,5;
9 - 9,6; хг = і3; Р = і
Схема «крест» условно устойчива, если соблюдено условие Ку-ранта-Фридрихса-Леви:
Ах <Л
г = — < і.
А5
Неявная схема
-і- (еп+і - 2еп + еп-і ї (еп+і - 2еп+і + еп+і ї
АхД Єі 2Єі + Єі ) А^Д Єі+і 2Єі +
(29)
Неявная схема с весами
і / п + і ГЬ п . п-і \ і / п + і Ґл \ п , п-і \
Х2(е - 2еі + е ) = А^2(стеі+і -(і - 2а)е +аеі-і ).
Ах
В двух последних случаях решаем разностную задачу методом прогонки.
Жесткое условие вида
Ф(т) =
I1, т<тс, ¡0, т>І0
(30)
особенно четко проявляет достоинства и недостатки разностных схем. Точное аналитическое решение в этом случае имеет вид
е, = А(Ж (т-5,)-Ж (Т-Т0-5,)),
(31)
[0, д < 0, д > 0,
Если г < г* = 1, как и в случае телеграфного уравнения, в местах разрыва наблюдаются нефизические осцилляции (рис. 4, а), которые практически исчезают, если г ^ г* (рис. 5, а).
а б
Рис. 4. Сравнение численного (квадратики) и аналитического (сплошная линия) решений волнового уравнения с прямоугольным (а) и синусоидальным (б) импульсом. Явная схема «крест». Момент времени Т : 0,625 • 10-2; г = 0,9. Время действия импульса т0 = 0,004
Если импульс имеет синусоидальную форму
Ф(т) =
Sin
0,
с \ п
—т
Чт0 J
, т<т0
т>т0,
(32)
осцилляций не наблюдается; численное решение (квадратики на рис. 4, б) хорошо согласуется с аналитическим решением (сплошная линия) для любого г < г* (рис. 4 и 5, б).
а б
Рис. 5. Сравнение численного (квадратики) и аналитического (сплошная линия) решений волнового уравнения с прямоугольным (а) и синусоидальным (б) импульсом. Явная схема «крест». Момент времени Т : 0,625 -10~2, г = 1,0
Неявная схема приводит к размытию решения в любом случае (рис. 6, а, б).
а б
Рис. 6. Сравнение численного (квадратики) и аналитического (сплошная линия) решений волнового уравнения с прямоугольным (а) и синусоидальным (б) импульсом. Неявная разностная схема. Момент времени Т : 0,625 -10_2
Заметим, что в случае неявной схемы наблюдается асимметрия численного решения относительно аналитического (передний фронт сглажен сильнее заднего). Это связано с тем, что время участия в численном расчете области переднего фронта больше, а следовательно, больше и накапливаемая погрешность схемы. С течением времени эта особенность устраняется, и профиль численного решения выравнивается относительно точного решения.
Явная разностная схема позволяет получить результаты, в большей степени согласующиеся с аналитическими. Но стоит отметить, что схема «крест» при варьировании временного шага дает значительное количество осцилляций в местах разрыва, что нежелательно и не имеет
физического смысла. Также необходимо, чтобы было соблюдено условие Куранта.
Таким образом, для решения связанной задачи была выбрана явная схема «крест».
5. Результаты численного исследования связанной модели
Как и в случае простых вариантов, были исследованы задачи двух типов.
Первая - с граничным условием для концентрации первого рода (соответствует, например, насыщению поверхности материала частицами из окружающей плазмы). В момент прекращения действия импульса длительностью т0 ставится условие отсутствия источников и
стоков массы на свободной поверхности: I = 0.
Вторая - с граничным условием для потока второго рода (например, обработка поверхности материала потоком частиц).
Для второго уравнения в любом случае на границе ставится условие равенства напряжений нулю, что приводит к соотношению
5 = 0: е = у(С - Со).
Решение связанной задачи для случая, если на левой границе задана концентрация, представлено на рис. 7.
Скорость распространения обеих волн - концентрационной и упругой - конечна. Но так как \ >> ^, упругая волна убегает вперед. Положение переднего фронта волны хорошо видно на рис. 7, б. Смене знака деформаций соответствует координата переднего фронта диффузионной волны. Расстояние между волнами с течением времени увеличивается. Уже к моменту времени т = 20 распределение концентраций вблизи свободной поверхности практически не изменяется (рис. 7, в), в то время как передний фронт диффузионной волны все же продвигается вперед. На рисунке показана только область вблизи свободной границы, где кривая деформации подобна кривой концентрации (рис. 7, г).
Решение для этих моментов времени можно назвать квазистати-ческим, если не принимать во внимание, что качественные особенности волн, отмеченные для рис. 7, а, б, имеют место вне области, изображенной на рис. 7, в, г, при 5 > 5,5 и 5 > 20.
0 1 2 3 4 £
в г
Рис. 7. Пример решения связанной задачи. Прямоугольный импульс. Моменты времени (а, б) т: 1 - 0,8, 2 - 4,0, 3 - 6,0; (в, г) т > 20 . Параметры модели ю = 100; у = -0,05; тг = 13,0. Явная схема. Время действия импульса т0 = 200тг
Пример решения задачи для прямоугольного импульса конечной длительности показан на рис. 8. Численное решение хорошо передает скорость распространения волн, положение фронтов и влияние волн друг на друга. Осцилляции уменьшаются с уменьшением шагов, так что физической природы они не имеют. Рисунки слева соответствуют моментам времени, меньшим времени импульса. Кривым с номерами 5, 6 (справа) - соответствуют моменты времени т > Т0. Особенности волн, отмеченные выше, сохраняются.
Пример решения задачи с граничными условиями второго рода представлен на рис. 9 для синусоидального импульса, т0 = 200тг. К моменту времени т упругая волна, как и в предыдущей задаче, доходит до точки £, = т . Деформации значительно меньше. Координата смены знака деформациями соответствует переднему фронту концентрационной волны.
Для т0 =тг качественная картина изменяется (рис. 10).
н-----■---1--------1----■---1---------1---■----1--------1---'---- . I ■ I ■ I ' I ■ . ■ I ' I ' I ■ I ' I „
О 1 2 3 4 5 6 £ О 2 4 6 8 10 12 14 16 18 20 5
б
Рис. 8. Пример решения связанной задачи. Прямоугольный импульс. Моменты времени Т : 1 - 0,8; 2 - 4,0; 3 - 6,0; 4 - 8,4; 5 - 15,0; 6 - 16,0; параметры модели: ю = 100; у = -0,05; тг = 13,0 . Явная схема. Время действия импульса т0 = 10,0
0,0 0,5 1,0 1,5 § 0 2 4 £
а
—0,0004-■—I—■—I—»—I—т—|—.—|—.—|—т— 0 0012
о 1 2 3 4 5 6 £, ’ 0 2 4 6 8 10 12 14 16 18 20 £,
б
Рис. 9. Пример решения связанной задачи. Синусоидальный импульс. Моменты времени Т : 1 - 0,8; 2 - 4,0; 3 - 6,0; 4 - 8,4; 5 - 15,0; 6 - 16,0. Параметры модели: ю = 100; у = -0,05; тг = 13,0 . Явная схема. Время действия импульса т0 = 200тг
Рис. 10. Пример решения связанной задачи. Синусоидальный импульс. Моменты времени Т : 1 - 1,25; 2 - 2,5; 3 - 3,75; 4 - 5,0; 5 - 6,25; 6 - 7,5; 7 - 10,0; 8 - 12,5; 9 -15,0. Параметры модели ю = 100; у = -0,05; тг = 13,0. Явная схема. Время действия
импульса т0 = 13,0
Из слабовогнутого фронт становится выпуклым. За постепенным «вхождением» волны через свободную поверхность следует ее отрыв от поверхности и распространение в глубь слоя.
Расчеты показывают, что все физические параметры оказывают влияние на форму волн и характер их эволюции во времени. Но всякий раз при переходе от одной области изменения параметров к другой требуется тщательный выбор параметров разностной схемы.
Заключение
Таким образом, в работе предложена изотермическая динамическая модель начальной стадии процесса внедрения ионов в поверхностный слой материала в условиях поверхностной обработки. Модель учитывает диффузию ионов, конечность времени релаксации потока массы, появление напряжений вследствие изменения состава поверхностного слоя и явление переноса компонента механической волной. Описан выбор безразмерных переменных и физический смысл получаемых при этом безразмерных комплексов. На основе анализа известных частных уравнений и классических разностных методов их чис-
ленного решения выбран метод численной реализации связанной модели. Выявлены качественные особенности концентрационной и упругой волн, взаимодействующих между собой. Смена качественного характера распределения деформаций в упругой волне связана с конечностью скорости распространения диффузионной волны.
Работа выполнена при финансовой поддержке государственного контракта 16.740.11.0122 в рамках федеральной целевой программы «Научные и научно-педагогические кадры инновационной России на 2009-2013 годы» по лоту «Проведение научных исследований научными группами под руководством докторов наук по следующим областям: математика; механика; информатика», шифр «2010-1.2.1-102-017».
Библиографический список
1. Ионная имплантация: сб. ст. / под ред. Дж.К. Хирвонена. - М.: Металлургия, 1985. - 392 с.
2. Амирханов И.В. Численный расчет температурных эффектов в материалах при облучении их тяжелыми ионами высоких энергий в рамках уравнений теплопроводности для электронов решетки // Препринт ОИЯИ. - Дубна, 2004. - С. 23. - ИКЬ: http://wwwinfo.jinr.ru/~ze-те1^гапЮ1/тёех.Мт1.
3. Амирханов И.В. Численное моделирование испарения металлов под действием импульсных ионных пучков // Препринт ОИЯИ. -Дубна, 2003. - С. 21. - ИЯЬ: http://wwwinfo.jinr.ru/~zeme1/grant01/in-dex.html.
4. Численное моделирование фазовых переходов в металлах, облучаемых импульсными пучками ионов // Препринт ОИЯИ / И.В. Амирханов, Е.В. Земляная, И.В. Пузынин [и др.]. - Дубна, 2001. - С. 15. -иКЬ: http://jdsweb.jinr.ru/record/57972/fi1es/p11-2001-164.pdf.
5. Бойко В.И., Евстигнеев В.В. Введение в физику взаимодействия сильноточных пучков заряженных частиц с веществом. - М.: Энер-гоатомиздат, 1988. - 136 с.
6. Взаимодействие импульсных пучков заряженных частиц с веществом / В.И. Бойко, В.А. Скворцов, В.Е. Фортов, И.В. Шаманин. -М.: Физматлит, 2003. - 288 с.
7. Модифицирование и легирование поверхности лазерными, ионными и электронными пучками / под ред. Дж. М. Поута, Г. Фоти,
Д.К. Джекобсона, А.А. Углова / пер. с англ. Н.К. Мышкина, А.В. Белого, В.М. Анищика. - М.: Машиностроение, 1987. - 424 с.
8. Князева А.Г. Нелинейные модели деформируемых сред с диффузией // Физическая мезомеханика. - 2011. - № 6. - 35-51 с.
9. Титановые сплавы в машиностроении / под ред. Г.И. Капыри-на. - Л.: Машиностроение, 1977. - 125 с.
10. Алюминий: свойства и физическое металловедение / под ред. Дж.Е. Хэтча. - М.: Металлургия, 1989. - 425 с.
11. Лыков А.В. Тепломассообмен: справочник. - М.: Энергия, 1978. - 450 с
12. Калиткин Н.Н. Численные методы. - М.: Наука, 1978. - 512 с.
References
1. Ion implantation / Ed. Hirvonen G.K. Treatise on Material Science and Technology. London: Acad. Press, 1980, 502 p.
2. Amirhanov I.V. Chislennij raschet temperaturnih effectov v materi-alah pri obluchenii ih tayzhelimi ionami visokih energij v ramkah uravnenij teploprovodnosti dlay electronov reshetki [The numerical calculation of thermal effects in materials irradiated by high-energy heavy ions in the heat equations for electrons lattice]. Preprint OIYAI, Dubna, 2004, pp. 23, available at: http://wwwinfo.jinr.ru/~zemel/grant01/index.html.
3. Amirhanov I.V. Chislennoe modelirovanie ispareniya metallov pod dejstviem impul’snih ionnih puchkov [Numerical simulation of the evaporation of metals under the action of pulsed ion beams]. Preprint OIYAI, Dubna, 2004, pp. 21, available at: http://wwwinfo.jinr.ru/~zemel/grant01/index.html.
4. Amirhanov I.V. Chislennoe modelirovanie fazovih perehodov v metal-lax, obluchaeih impul’snimi puchkami ionov [Numerical modeling of phase transitions in metals irradiated by pulsed ion beams]. Preprint OIYAI, Dubna, 2001, pp. 15, available at: http://jdsweb.jinr.ru/record/57972/files/p11-2001-164.pdf.
5. Boyko V.I., Evstigneev V.V. Vvedenie v fiziku vzaimodejstviay sil’notochnih puchkov zarayzhennih chastic s vewestvom [Introduction to the physics of the interaction of intense beams of charged particles with matter]. Мoscow: Energoatomisdast, 1988, 136 p.
6. Boyko V.I., Skvorcov V.A., Fortov V.E., Shamanin I.V. Vzai-modejstvie impul’snih puchkov zarayzhennih chastic s vewestvom [Interaction of charged particle beams with matter]. Мoscow: Fizmatlit, 2003, 288 p.
7. Surface modification and alloying by laser, ion, and electron beams. Ed. J. M. Poate, G. Foti, D. C. Jacobson. NATO Scientific Affairs Division, Plenum Press, 1987, 424 p.
8. Knayzeva A.G. Nelineinie modeli deformiruemih sred s diffuziey [Nonlinear models of strained media with diffusion]. Fizicheskaya mezome-hanika, 2011, no. 6, рp. 35-51.
9. Titanovie splavi v mashinostroenii [Titanium alloys in engineering]. Ed. G.I. Kapirina. Leningrad: Mashinostoenie, 1977, 125 p.
10. Aluminij: svojstva i fizicheskoe metallovedenie [Aluminum: properties and physical metallurgy]. Ed. D.E. Hatch. Мoscow: Metallurgiay, 1989,425 p.
11. Likov A.V. Teplomassoobmen [Heat and Mass Transfer]. Мoscow: Energiay, 1978, 450 p.
12. Kalitkin N.N.Chislennie metodi [Numerical methods]. Мoscow: Nauka, 1978, 512 p.
Об авторах
Ильина Елена Сергеевна (Томск, Россия) - инженер лаборатории моделирования физико-химических явлений в современных технологиях Национального исследовательского Томского Политехнического Университета (634050, г. Томск, пр. Ленина, 30, е-mail: Lina-sergg@mail.ru).
Демидов Валерий Николаевич (Томск, Россия) - кандидат физико-математических наук, доцент Кафедры физики высоких технологий в машиностроении Национального исследовательского Томского Политехнического Университета (634050, г. Томск, пр. Ленина, 30, е-mail: vn_demidov@mail .ru).
Князева Анна Георгиевна (Томск, Россия) - доктор физ.-мат. наук, профессор, зав. лабораторией моделирования физико-химических явлений в современных технологиях Национального исследовательского Томского Политехнического Университета (634050, г. Томск, пр. Ленина, 30), главный научный сотрудник Института физики прочности и материаловедения СО РАН (634021, Томск, Академический пр. 2/4, e-mail: anna-knyazeva@mail.ru).
About the authors
Il’ina Elena Sergeevna (Tomsk, Russian Federation) - engineer of laboratory of Modeling of physical-chemical processes in modern technologies, National Research Tomsk Polytechnic University (30, Lenina av., 634050, Tomsk, Russian Federation, e-mail: Linasergg@mail.ru).
Demidov Valeriy Nikolaevich (Tomsk, Russian Federation) - Ph. D. of Physical and Mathematical Sciences, Ass. Professor, Department of High Technology in Engineering Industry, National Research Tomsk Polytechnic University (30, Lenina av., 634050, Tomsk, Russian Federation, e-mail: vn_demidov@mail.ru).
Knyazeva Anna Georgievna (Tomsk, Russian Federation) - Doctor of Physical and Mathematical Sciences, Professor, head of laboratory of Modeling of physical-chemical processes in modern technologies National Research Tomsk Polytechnic University (30, Lenina av., 634050, Tomsk, Russian Federation,); principal researcher of Institute of Strength Physics and Material Science (2/4, Academicheskii av., 634021, Tomsk, Russian Federation, e-mail: anna-knyazeva@mail.ru).
Получено 14.09.2012
