Методические указания к выполнению домашнего задания по курсу «Методы моделирования в материаловедении»

Московский государственный технический университет 
имени Н.Э. Баумана 

В.И. Третьяков, А.Ю. Ампилогов, М.А. Хасянов 

Методические указания  
к выполнению домашнего задания 
 по курсу «Методы моделирования  
в материаловедении» 

М о с к в а 
Издательство МГТУ им. Н.Э. Баумана 
2 0 0 6 

УДК 620.22(076.5) 
ББК 30.3 
Т66 

Рецензент В.Н. Симонов 

Третьяков В.И., Ампилогов А.Ю., Хасянов М.А. 
Т66 
 
Методические указания к выполнению домашнего задания 
по курсу «Методы моделирования в материаловедении». – М.: 
Изд-во МГТУ им. Н.Э. Баумана, 2006. – 59 с.: ил.  

ISBN 5-7038-2858-9 

Предложены методология системного анализа, методы физического 
моделирования для решения задач моделирования в машиностроительном 
материаловедении. Рассмотрены общие алгоритмы 
создания математических моделей, планирования и проведения 
численного эксперимента с анализом результатов моделирования. 
Приведены примеры моделирования процессов диффузионного насыщения 
при химико-термической обработке и тепловых процессов 
нагрева-охлаждения при термической обработке. 
Для студентов специальности «Материаловедение» (МТ-8), 
изучающих курсы «Основы автоматизированного проектирования» 
и «Методы моделирования в материаловедении» в МГТУ 
им. Н.Э. Баумана. 
Ил. 34. Табл. 8. Библиогр. 6 назв. 

                УДК 620.22(076.5) 
                 ББК 30.3 

ISBN 5-7038-2858-9  
 МГТУ им. Н.Э. Баумана, 2006 

ВВЕДЕНИЕ 

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

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

1. ТЕОРЕТИЧЕСКАЯ ЧАСТЬ 

1.1. Математическое моделирование диффузионных процессов 

при химико-термической обработке 

Уравнения диффузии 

Моделирование диффузионных процессов в металле основывается 
на уравнениях Фика: 
первом 

 
c
J
D x

∂
=−
∂  
(1) 

и втором  

 

2

2 .
c
c
D

x

∂
∂
=
∂τ
∂
 
(2) 

Здесь D – коэффициент диффузии; c – концентрация; τ – время; x – 
координата. 
Первое уравнение описывает удельный поток диффундирующего 
элемента в металле. Знак «минус» означает, что поток направлен 
из области с большей концентрацией в область с меньшей 
концентрацией.  
Второе уравнение Фика описывает изменение концентрации 
диффундирующего вещества c(x, τ) в пространстве и во времени. 
Это уравнение непосредственно следует из уравнения баланса 
вещества при диффузии и выражения для потока. Оба 
уравнения описывают одномерную диффузию вдоль оси x; в 
качестве допущения для упрощения решения поставленной задачи 
принимаем, что коэффициент диффузии D не зависит от 
концентрации диффундирующего элемента, а только от температуры 
процесса. 

Начальные и граничные условия 

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

 
0
( , 0)
( ),
c x
c
x
=
 
(3) 

а в случае постоянной начальной концентрации – условием 

 
0
( ,0)
const.
c x
c
=
=
 
(4) 

Если рост диффузионного слоя контролируется диффузией в 
металле, то используется граничное условие 1-го рода 

 
(0, )
( ),
S
c
c
τ =
τ  
(5) 

т. е. концентрация на поверхности детали cS (при x = 0) является 
функцией времени или, в частном случае, постоянна: 

 
(0, )
const.
S
c
c
τ =
=
 
(6) 

Концентрация на поверхности детали определяется, как правило, 
экспериментальным путем. 
Когда рост диффузионного слоя контролируется процессами, 
происходящими на поверхности металла, например адсорбцией 
или химической реакцией, используется граничное условие 2-го 
рода. В этом случае при решении уравнения (2) задается поток вещества 
через поверхность металла как функция времени. 
Граничное условие 3-го рода наиболее достоверно описывает 
реальный процесс массопереноса на поверхности металла при химико-
термической обработке: 

 
(
) 0
S
c
D
k c c
x

∂
−
+
−
=
∂
 для x = 0, 
(7) 

где k – константа скорости процесса, происходящего на поверхности 
детали. 

Аналитическое решение уравнения диффузии 

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

2

2
С
C
D

x

∂
∂
=
∂τ
∂
 

и начальными и граничными условиями, имеющими вид 

τ = 0, С = С1 для всех x > 0, 

τ = 0, С = С2 для всех x = 0, 

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

( )
.
x
С
t
λ
=
 

Тогда 

3 2
λ
,
λ
λ 2
C
dC
dC
x
t
d
t
d
t
∂
∂
=
=−
∂
∂
 

1 ,
С
C
C
x
x
t

∂
∂
∂λ
∂
=
=
∂
∂λ
∂
∂λ
 

2
2

2
2
1
1.
C
C
C

x
t
x
t

∂
∂
∂
∂λ
∂
 
 
=
=
 
 
 
 
∂λ
∂λ
∂
∂
∂λ
 

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

2

3 2
2
1.
2
С
x
C
D
t
t
∂
∂
−
=
∂λ
∂λ
 

Заметим, что 
3 2
1 1.
2
2
x
t
t
=
λ
 После подстановки этого выражения 

в последнее дифференциальное уравнение получим 

2

2
2
.
dC
d C
D
d
d
λ =−
λ
λ
 

Первую производную dC/dλ представим в виде 

.

n
a
dC
Ae
d

− λ
=
λ
 

Подставив это выражение в предыдущее дифференциальное 
уравнение, получим 

1
2
.

n
n
a
a
n
Ae
DAe
an
− λ
− λ
−
λ=
λ
 

Приравняв значения показателей степени экспоненты с левой и 
правой сторон уравнения, а также показатели степени при переменной 
λ, получим n = 2, а = 1/(4D). Подставив эти значения в 
предполагаемое решение для первой производной, получим 

2 /(4
),
D
dc
Ae
d

−λ
=
λ
 

где А – постоянная интегрирования. 
Повторное интегрирование позволяет получить общее решение 
в виде 

2 /(4
)

0
.
D
C
A e
d
B
λ −λ
=
λ+
∫
 

Для удобства интегрирования можно ввести новую переменную  

/(2
).
D
ξ=λ
 

Тогда общее решение будет иметь следующий вид: 

2
2
/(2
)
/(2
)
'"

0
0
/(2
)
.
D
x
Dt
C
A
D
e
d
B
A
e
d
B
λ
−ξ
−ξ
=
ξ+
=
ξ+
∫
∫
 

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

2

0
.
2
e
d
∞ −ξ
π
ξ=
∫
 

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

при t = 0    C = C1 для 
,
x→∞  

при t = 0     C = C2 для 
0.
x=
 

Подставляя условия, определяющие граничные и начальные 
значения решения в общее уравнение, определим постоянные интегрирования: 


1
,
2
C
A
B
π
=
+
 

С2 = В. 

Тогда  

2
1
(
)2,
C
C
A
−
=−
π
 

В = С2, 

2
/(2
)

2
2
1
0

2
(
)
.

x
Dt
C
C
C
C
e
d
−ξ
=
−
−
ξ
∫
π
 

Таким образом, при диффузии в полуограниченный образец  
(0 < x < ∞, где x – расстояние от поверхности) с нулевой начальной 
концентрацией c (x, 0) = c0 = 0 через поверхность (x = 0), на которой 
поддерживается постоянная концентрация cS, не зависящая от времени 
c (0, τ) = cS, распределение концентрации описывается уравнением 

( , )
1 erf
,
2
S
x
c x
c
D
⎛
⎞
τ =
−
⎜
⎟
⎝
⎠
τ
 

где функция ошибок erf определяется как 

2

0

2
erf
exp(
)
.
z

z
d
=
− ξ
ξ
∫
π
 

Решение диффузионной задачи численным методом 
 конечных разностей 

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

Рис. 1. Пространственно-временная сетка 

В пространственной области выбирают узлы пространственной 
сетки: конечное число значений координат x1, x2, ..., xN (N – максимальное 
значение индекса вдоль оси координат); для временной 
переменной выбирают конечное число значений τ0, τ1, ..., τJ (узлы 
временной сетки; J – максимальное значение индекса вдоль оси 
времени). Узлы пространственной сетки располагаются на равном 
расстоянии (рис. 2): 

1
n
n
x
x
x
+ −
=∆  для n = 1, N. 

Рис. 2. Пространственная сетка 

На втором этапе для каждого узла пространственной сетки записывают 
уравнение (2) относительно значений искомой функции в узловых 
точках. Члены, входящие в уравнения (2), аппроксимируют, выражая 
их через значения сеточной функции cn,j, являющейся лишь приближенным 
значением функции концентрации c (x, τ). Пренебрегая 
разностью между значением сеточной функции cn,j и функции концентрации 
c (x, τ), для временной производной используют выражение 

 
,
1
,

,
;

n
j

n j
n j

x

c
c
c
+

τ

−
∂
=
∂τ
∆τ
 
(8) 

для первой производной по координате x – выражение 

 
1,
,

,
;

n
j

n
j
n j

x

c
c
c
x
x

+

τ

−
∂
=
∂
∆
 
(9) 

для второй производной по координате x – выражение 

 

2
1,
,
1,
2
2
,

2
.

n
j

n
j
n j
n
j

x

c
c
c
c
x
x

+
−

τ

−
+
∂
=
∂
∆
 
(10)