Компьютерное моделирование нанотехнологий, наноматериалов и наноструктур : диффузия

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РФ 

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ  
ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ  
«НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ ТЕХНОЛОГИЧЕСКИЙ УНИВЕРСИТЕТ «МИСиС» 

 

 
 
 

 

 

 

 
 

 

№ 1894 

Кафедра полупроводниковой электроники и физики полупроводников
 
Ю.В. Осипов 
М.Б. Славин 
 

Компьютерное моделирование
нанотехнологий, наноматериалов
и наноструктур 

Диффузия 

Учебное пособие 

Рекомендовано редакционно-издательским 
советом университета 

Москва 2011 

УДК 621.315.592 
 
О-74 

Р е ц е н з е н т  
д-р техн. наук, проф. Г.Д. Кузнецов 

Осипов, Ю.В. 
О-74  
Компьютерное моделирование нанотехнологий, наноматериалов и наноструктур : диффузия : учеб. пособие / Ю.В. Осипов, 
М.Б. Славин. – М. : Изд. Дом МИСиС, 2011. – 73 с. 
ISBN 978-5-87623-420-9 

В учебном пособии представлены выводы уравнения диффузии и уравнения Пуассона. Получено уравнение диффузии, учитывающее наличие электрического поля в области контакта двух различных материалов. Рассмотрено влияние химических реакций на процессы диффузии. Приведены решения 
диффузионных задач, в которых учтено: отличие коэффициентов диффузии 
на различных границах раздела фаз, присутствие внутреннего электрического поля, наличие химических реакций в области контакта материалов. Представлены программы на языке С++, дающие возможность читателю самостоятельно решать задачи, указанные выше, и их комбинации. 
Соответствует государственному образовательному стандарту дисциплины «Компьютерное моделирование нанотехнологий, наноматериалов и наноструктур». 
Предназначено для магистров, обучающихся по направлениям 210100 
«Электроника и наноэлектроника» и 150100 «Материаловедение и технологии материалов». 
УДК 621.315.592 

ISBN 978-5-87623-420-9 
© Осипов Ю.В., 
Славин М.Б., 2011 

ОГЛАВЛЕНИЕ 

Введение................................................................................................4 
1. Уравнение диффузии........................................................................5 
2. Решение уравнения диффузии ......................................................11 
3. Уравнение Пуассона.......................................................................31 
4. Уравнения диффузии при наличии электрического поля...........41 
5. Решение уравнения диффузии при наличии электрического  
    поля ..................................................................................................45 
6. Влияние химических реакций на процессы диффузии...............61 
Библиографический список...............................................................72 

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

Фейнман Р., Лейтон Р., Сэндс М.  
Фейнмановские лекции по физике 

Введение 

В термодинамически равновесных системах, как известно, интенсивные параметры не зависят от пространственных координат,  
например 

 
grad T = 0, grad P =0, 
(1) 

где T – температура;  
P – давление.  

Если условие (1) не выполняется, то в системе возникают необратимые процессы переноса. В случае появления в системе потоков 
тепла говорят о теплопроводности, при переносе зарядов – об электропроводности. Явление, связанное с перемещением вещества 
вследствие теплового движения, называют диффузией (от лат. 
diffusio – распространение, растекание). Необходимо обратить внимание читателя на очевидное, но иногда игнорируемое обстоятельство принадлежности понятия «диффузия» к науке, которая называется 
«Неравновесная термодинамика». Не считаться с этим фактом нельзя, потому что это может привести к неправомерным выводам при 
обсуждении результатов экспериментального исследования или теоретического изучения рассматриваемого явления, например, к появлению «равновесного» коэффициента самодиффузии [1]. Необходимо подчеркнуть, что для исследования любой проблемы, решения 
любой задачи, прежде всего, следует обусловить сам объект исследования, то есть определить пространство, в котором будут рассмотрены проблемы, методы и способы работы с параметрами – координатами пространства. В данном пособии диффузия будет рассматриваться в рамках линейной неравновесной термодинамики. Это означает, что, хотя в целом состояние системы неравновесно, отдельные 
ее малые части квазиравновесны и имеют медленно изменяющиеся 
во времени и от точки к точке интенсивные термодинамические параметры. Размеры этих малых равновесных частей неравновесной 
системы и времена изменения термодинамических параметров в них 

определяются экспериментально. Свойства неравновесной системы 
при этом определяются локальными термодинамическими параметрами, которые зависят от пространственных координат и времени 
только через характеристические термодинамические параметры, для 
которых справедливы уравнения термодинамики. Причем при малых 
отклонениях от равновесия связь между потоками Ii и термодинамическими силами Xk (интенсивными параметрами термодинамической 
системы, такими как температура, давление, напряженность электрического поля, химический потенциал и т.д.) 1 линейна: 

 

1

n

i
ik
k

k
L

=
=∑
I
X . 
(2) 

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

1. Уравнение диффузии 

Рассмотрим систему, содержащую k компонентов, при этом допустим, что i-й компонент имеет градиент химического потенциала, 
который, в свою очередь, обусловлен градиентом концентрации этого компонента. По определению поток1 ji(x,y,z) равен количеству частиц, которые пересекают в единицу времени поверхность единичной 
площади, то есть 

 
( , , , )
( , , , )
( , , , )
i
i
i
x y z t
c x y z t
x y z t
=
j
v
,  
(3) 

где x, y, z 
– 
пространственные координаты;  
t 
 
– 
временная координата;  
ci(x, y, z) 
– 
концентрация i-го компонента;  
vi(x, y, z) 
– 
средняя скорость частиц, которую выразим через  
  
 
 
подвижность ui(x, y, z) и силовую характеристику, 
  
 
 
вызвавшую поток, в системе fi(x, y, z) следующим 
  
 
 
образом: 

––––––––– 
1 Векторные величины в пособии выделены жирным шрифтом. 

vi (x, y, z, t) = 
)
,
,
,
(
)
,
,
,
(
t
z
y
x
t
z
y
x
u
i
i
f
.  
(4) 

Прежде чем продолжить вывод уравнения диффузии, необходимо 
сделать небольшое отступление, связанное с физическим смыслом 
величины fi(x,y,z). 
Из курса общей физики известно, что если в пространстве определено распределение электрического потенциала φ(x, y, z), то в этом 
пространстве определена и его силовая характеристика Ei(x, y, z) – 
напряженность электрического поля, причем 

 
( , , )
i x y z
E
= – grad ( , , )
x y z
ϕ
  
(5) 

или 

 
( , , )
i x y z
E
= –∇
( , , )
x y z
ϕ
, 
(5a) 

где ∇  – оператор набла, который в декартовой системе координат 
        имеет вид 

 
∇ =
x
y
z

∂
∂
∂
+
+
∂
∂
∂
i
j
k , 
(5б) 

где, в свою очередь, i, j, k – ортогональные единичные векторы, определяющие декартову систему координат.  

Справедливость равенства правых частей выражений (5) и (5а) 
легко показать. Действительно, 

 
∇
( , , )
x y z
ϕ
=
x
y
z

⎛
⎞
∂
∂
∂
+
+
⎜
⎟
∂
∂
∂
⎝
⎠

i
j
k
( , , )
x y z
ϕ
=

( , , )
( , , )
( , , )
x y z
x y z
x y z
x
y
z
∂ϕ
∂ϕ
∂ϕ
=
+
+
=
∂
∂
∂
i
j
k
grad
( , , )
x y z
ϕ
. 

Здесь же следует напомнить, что скалярное произведение оператора 
набла и некоторого вектора b(x,y,z) = bx(x,y,z)i + by(x,y,z)j + bz(x,y,z)k 
определяет скалярную величину, которую называют дивергенцией: 

∇
=
)
,
,
(
z
y
x
b
(
)
( , , )
( , , )
( , , )
x
y
z
b x y z
b x y z
b x y z
x
y
z
⎛
⎞
∂
∂
∂
+
+
+
+
=
⎜
⎟
∂
∂
∂
⎝
⎠

i
j
k
i
j
k

( , , )
( , , )
( , , )
y
x
z
b
x y z
b
x y z
b x y z
x
y
z

∂
⎛
⎞
∂
∂
=
+
+
=
⎜
⎟
∂
∂
∂
⎝
⎠

div ( , , )
x y z
b
. 

Обратим внимание на следующие словосочетания: электрический 
потенциал, термодинамический потенциал, потенциал Гиббса. Общим для этих выражений является слово «потенциал». Этот факт не 
случаен: подвергая указанные выше параметры однотипным преобразованиям, мы должны получать однотипные по физическому 
смыслу результаты (величины), то есть, воздействуя оператором набла на электрический потенциал, мы получаем, как было показано 
выше, характеристику поля, являющуюся движущей силой процессов, происходящих с заряженными частицами. Аналогично, действуя 
оператором набла на термодинамический потенциал, мы непременно 
должны получить подобный параметр. При этом надо помнить, что в 
термодинамике количества теплоты и работы являются функциями 
процессов (адиабатного, изохорического, изобарического и т.д.). Поэтому, например, в изолированной системе при квазиизохорическом 
процессе за все изменения ответственен градиент температуры: 
∇(dU) = ∇(∂U/∂S)V,N dS = ∇TdS, а в квазиадиабатном – градиент давления: ∇(dU) = ∇(∂U/∂V)S,N dV = –∇PdV (здесь U – внутренняя энергия; 
S – энтропия; V – объем; N – число частиц в системе; P – давление). 
В нашем случае, как мы условились ранее, поток i-го компонента 
обусловлен градиентом химического потенциала μi, поэтому1 

 
if  = – ∇
iμ .  
(6) 

Подставив (6) в (4), а результат подстановки в (3), получим 

 
ij  = – 
ic
iu ∇
iμ .  
(7) 

Собственно говоря, (7) и есть первый закон Фика, записанный в 
самом общем виде и справедливый как для реальной, так и для идеальной систем. А теперь установим зависимость изменения количества i-го компонента Ni от времени t в некотором локальном объеме 
∆V в результате истечения частиц через поверхность ∆S, ограничивающую ∆V: 

 
d
i
i
S

N
t
Δ

∂
= −
∂
∫ j S . 
(8) 

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

В уравнении (8) знак минус поставлен вследствие того, что количество частиц в рассматриваемом объеме уменьшается со временем 
и, соответственно, левая часть последнего выражения должна быть 
меньше нуля, а под dS имеется в виду 

 
d
dS
=
S
n , 
(9) 

где n  – единичный вектор, направленный нормально к элементу поверхности dS изнутри исследуемой области.  

Зная пространственное распределение компонентов системы 
ci(x,y,z), легко рассчитать Ni: 

 
d
i
i
V

N
c V

Δ
= ∫
.  
(10) 

Подставим в (10) в (8), получим 

 
d
i
V

c V
t Δ

∂
∂ ∫
d
i
S
Δ
= −∫ j S .  
(11) 

Используя теорему Остроградского – Гаусса, в правой части (11) 
перейдем от интегрирования по поверхности к интегрированию по 
объему: 

 
d
i
V

c V
t Δ

∂
∂ ∫
div
d
i
V

V

Δ
= − ∫
j
.  
(12) 

Поскольку область интегрирования в (12) выбрана произвольно и 
ее размер в левой и правой частях одинаков, равенство будет выполняться в случае, когда 

 
ic
t

∂
∂
div i
= −
j ,  
(13) 

или, используя оператор набла, 

 
ic
t

∂
∂

i
= −∇j . 
(14) 

Формулы (13) – (14) выражают второй закон Фика. 
Но уравнения (7) и (13) – (14) в литературе не связывают с именем 
немецкого ученого Фика, поскольку указанные формулы получены 
для любых (идеальных и неидеальных систем), а Фиком сформули

рован закон диффузии для идеальных систем, в соответствии с которым «количество соли, проходящее за известное время в направлении убывающей концентрации через некоторый элемент поверхности, пропорционально величине этого элемента, промежутку времени, величине убывания концентрации на месте нахождения элемента 
поверхности по направлению течения»[3]. 
Покажем допущения, которые необходимо сделать, чтобы от (7) и 
(14) перейти к традиционным уравнениям диффузии, названными законами Фика. Имея в виду, что зависимость химического потенциала 

iμ  от температуры Т и активности i-го компонента 
ia  имеет вид1 

 
0
Б ln
i
i
i
k T
a
μ = μ
+
, 
(15) 

преобразуем (7) и (14) к 

 
Б
i
i

i
i

i

c u
k T
a
a
= −
∇
j
 
(16) 

и 

 
ic
t

∂
∂
=
Б
(
)
i
i

i

i

c u
k T
a
a
∇
∇
 
(17) 

соответственно. Применяя правила действий с дифференциальными 
операторами [4], получаем из (17) 

 
ic
t

∂
∂
 = 
2

Б
Б
i
i
i
i

i
i

i
i

c u
c u
k T
a
k T
a
a
a
⎛
⎞
∇
∇
+
∇
⎜
⎟
⎝
⎠
. 
 (18) 

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

––––––––– 
1 Во всех последующих формулах посредством kБ обозначена постоянная 
Больцмана. 

и то для достаточно ограниченного числа начальных и граничных 
условий. В случае идеальной системы выражение (7) трансформируется в уравнение вида 

 
Б
i
i
i
k Tu
c
= −
∇
j
. 
 (19) 

Три первых сомножителя в правой части (19) объединяют общим 
названием «коэффициент диффузии», то есть 

 
Б
i
i
D
k Tu
=
. 
 (20) 

В литературе равенство (20) известно как формула Эйнштейна. С 
учетом (20) уравнение (19) принимает ту форму записи, которую и 
называют первым законом Фика: 

 
i
i
i
D
c
= −
∇
j
, 
 (21) 

а уравнение (18) при дополнительном условии, что коэффициент 
диффузии не зависит от пространственных координат, принимает 
вид, который общепринято называть вторым  законом Фика или 
уравнением диффузии 

 
ic
t

∂
∂
=
Б
i
i
k Tu
c
Δ
. 
 (22) 

С учетом уравнения (20) 

 
ic
t

∂
∂
=
i
i
D
c
Δ
. 
 (23) 

В выражениях (22) и (23) посредством символа ∆ (дельта) обозначен оператор Лапласа, который в декартовой системе координат имеет вид 

 

2
2
2

2
2
2
x
y
z

∂
∂
∂
Δ = ∇⋅∇ =
+
+
∂
∂
∂
. 
(24) 

Для того чтобы все же учесть неидеальность системы, зависимость коэффициента активности и подвижности от пространственных координат переносят на коэффициент диффузии. Схема получения уравнения диффузии в этом случае следующая. Уже в (7) полагают, что Di = Di(x,y,z) – функция пространственных координат x, y, z,  
а затем уравнение (7) подставляют в (17), которое принимает вид  

ic
t

∂
∂
=
(
)
i
i
D
c
∇
∇
. 
(25) 

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

2. Решение уравнения диффузии 

Рассмотрим следующую задачу. Дана пластина (рис. 1), на которой 
сформировали маскирующий слой (допустим, оксид кремния на кремниевой пластине), в котором, в свою очередь, образовали окно шириной 2l и длиной, соизмеримой с линейными размерами пластины. 
Диффузант поступает из источника неограниченной емкости. 
Преобразуем уравнение диффузии (23) для трехмерного декартова 
пространства – x, y, t, в случае идеальной однокомпонентной системы 

 
( , , )
c x y t

t
∂
∂
 = 

2
2

2
2

( , , )
( , , )
c x y t
c x y t
D
x
y

⎛
⎞
∂
∂
+
⎜
⎟
∂
∂
⎝
⎠
 
(26) 

X

Y
0

2 l

Область диффузанта

Область кристалла

Маска
Окно

 

Рис. 1. Схема процесса диффузии