CC BY
145
УДК 519.624.3
MPI реализация алгоритмов для 2D и 3D моделирования фазовых переходов в материалах, облучаемых тяжёлыми ионами, в рамках модели термического пика
И. В. Амирханов, Е. В. Земляная, Н. Р. Саркар, И. С. Сархадов, З. К. Тухлиев, З. А. Шарипов
Лаборатория информационных технологий Объединённый институт ядерных исследований ул. Жолио-Кюри, д.6, Дубна, Московская область, 141980, Россия
В работе представлена MPI реализация метода 2D и 3D расчётов эволюции температурных полей и динамики фазовых переходов, возникающих в материалах при облучении тяжёлыми ионами высоких энергий и импульсными ионными пучками. Используется модифицированная модель термического пика, основанная на системе двух связанных уравнений теплопроводности, описывающих тепловые процессы в электронной и ионной подсистемах облучаемого тяжёлыми ионами образца. Численное решение этих уравнений осуществляется в цилиндрической системе координат как в аксиально симметричном случае (2D), так и с учётом нарушения аксиальной симметрии в моделируемой системе (3D). Моделирование динамики фазовых переходов реализовано на основе задачи Стефана в рамках энтальпийного подхода. Представлена математическая постановка задачи, приведены формулы, определяющие конечно-разностную вычислительную схему; обсуждаются особенности параллельной компьютерной реализации на базе технологии MPI. Приведены численные результаты, подтверждающие эффективность разработанного подхода и соответствующей MPI/C+—+ программы. Показано, что результаты численного моделирования согласуются с известными экспериментальными оценками размеров треков, образующихся в облучаемых тяжёлыми ионами образцах.
Ключевые слова: моделирование, численные методы, фазовый переход, модель термического пика, параллельный алгоритм.
1. Введение
Одной из широко применяемых моделей для описания процессов, протекающих при взаимодействии заряженных частиц с материалами, является модель термического пика (МТП), базирующаяся на системе двух связанных уравнений теплопроводности для электронного газа и кристаллической решётки [1]. МТП использована для объяснения эффектов трекообразования в металлах, диэлектриках [2], полупроводниках [3] и в бесструктурных аморфных сплавах [4,5]. Результаты развития и практического применения данной модели можно найти в работах [6,7].
Нужно отметить, что на практике формально трёхмерная модель часто используется в одномерном приближении (например, в [5,8]). В работах [6,7] связанные уравнения теплопроводности в рамках МТП решаются сеточными методами в цилиндрической системе координат с учётом аксиальной симметрии (2D расчёты). Условие аксиальной симметрии выполняется для многих мишеней с большой точностью, однако существуют системы, для которых симметрия нарушается, что требует специального рассмотрения (3D случай) и приводит к существенному увеличению объёма вычислений. Это делает актуальной разработку параллельных алгоритмов и программ для проведения расчётов на многопроцессорных вычислительных системах.
Мы представляем параллельную реализацию алгоритма для численного решения уравнений МТП, модифицированной для исследования тепловых процессов с
Статья поступила в редакцию 31 мая 2013 г. Работа выполнена при финансовой поддержке Минобрнауки России в рамках государственного контракта №07.524.12.4019 от 17.05.2012.
учётом фазовых переходов в аксиально несимметричном случае. Работа построена следующим образом. В разделе 2 дана общая математическая постановка задачи в рамках МТП для несимметричной мишени с использованием энтальпийного метода моделирования фазовых переходов. В разделе 3 описывается вычислительная схема. В разделе 4 обсуждаются особенности параллельной программной реализации вычислительной схемы. Далее, в разделе 5 представлены результаты численного моделирования облучения несимметричной никелевой мишени ионами урана, демонстрирующие возможности разработанного подхода. Раздел 6 по-свящён облучению материалов импульсными пучками ионов. При всех отличиях с физической точки зрения данного процесса от облучения одиночными ионами, математическую постановку задачи можно свести к описанной в разделе 1. Существенно меняется лишь функция, «отвечающая» за форму источника облучения. Кроме того, в этом случае нами вводится нелинейная зависимость физических параметров мишени от температуры, что делает модель более реалистичной, хотя и усложняет численное исследование. Развитый подход применён для моделирования фазовых переходов при облучении образцов железа импульсными пучками углерода. В разделе 7 обсуждаются численные результаты в сравнении с известными экспериментальными оценками.
2. Учёт фазовых переходов в рамках модели термического
пика
В основе МТП лежит система двух связанных уравнений теплопроводности для электронного газа и кристаллической решётки. Эта система уравнений в цилиндрической системе координат имеет следующий вид:
дте (\ д ( дте\ 1 д( дтл
+§~z(Xe(Te))- g(Te)(Te -Tt) + Ae(r,^,z,t), (1)
Сг(Тг^ = (г* Ь^-Ж) + ^ду д^) +
+д (Хг(Тг) дТ)) +9(Те)(Те-Тг ) + Мг,<Р,*, V. (2)
Здесь ось г направлена вдоль пучка и перпендикулярна облучаемой поверхности; Те(г, у, г, ¿), Тг(г, у, г, 1) — температурные поля электронного газа и решётки облучаемого образца; Се, Сг
и Ае, Хг — соответственно, удельные теплоёмкости и коэффициенты теплопроводности электронного газа и решётки; — коэффициент электрон-фононного взаимодействия электронной подсистемы с решёткой; Ае(г, у, г, 1) и Аг(г, у, г, 1) — объёмные плотности энергии, вносимые ионом в электронную и решёточную подсистемы, соответственно.
Функции объёмной плотности энергии выбираются в факторизованном виде по всем аргументам, а именно:
Ае,г(г,у,z, г) = /1(0/2((р)h(t)Sineh■pЪ(z),
причём каждая функция /г, (г = 1, 2, 3) нормирована на единицу:
С 2" СО ^тах
I /1(г)г&г = 1,1 ¡2Ш<Р =1, ! = 1,1 (Sinel(z) + Spъ(z))dz = Е0,
где Е0 — энергия падающего иона, и 5рь соответственно, определя-
ют потери энергии иона на возбуждения электронов и на фононные возбуждения в решёточной подсистеме. Эти функции вычисляются программой 8ШМ-2012 (http://www.srim.org).
Функции fi, (г = 1,2,3) выбираем следующим образом:
fi(r) = -1 exp(-г /го), ¡2(р) = 21 ( 1+ cos тр + с2т sin ту)
/ м \
1 + (cim cos т<£ + С2т sin mip)
V т=1 J
f (Л 1 ( V - 5to)2
fs(t) = exp -
^ Ч 252
Здесь г о, to, St — физические параметры [9], М — число гармоник.
Предполагая, что физические параметры Ce,CJi,\e,\i, q постоянные, т.е. не зависят от Te,i (линейная задача), решение системы (1)—(2) будем искать в виде:
м
Te,i(r,<p,z,t) = To,e,i(r,z,t) + ^ Trn,e,í(r,z,t)(cim cos my + C2m sin m<p).
rn=1
Тогда для функции Тш,е,{(г, z, t), (m = 0,1, 2,..., M) получаем систему уравнений: dTm,,e
=
/1 д ( дТ \ т2Т д2Т \
= ЛЧ --^ + ^2г) - 9(Тт,- - + Am,e(r,z,t), (3)
С,
дТ„
dt
= д. (1 — (гдТтА
\ г дг \ дг )
т2Тг
д 2Т
+
dz2
+ g(Tm,e - Tm,i) + Am,i(r,Z,t). (4)
Система (3)—(4) решается независимо для каждого значения т со следующими начальными и граничными условиями:
To,e,í(r,z,t)lt=o = То, Tm,e,i(r,z,t)lt=o = 0, т = 1,2,...,М,
dTo,e,i(r,Z,t)
дг
= 0, Tm,e,i(r,z,t)lr=0 = 0, т = 1,2,.
г=0
dTm,e,i(r,z,t)
,М,
дг
dTm,e,i(r,z,t)
dz
0,
= 0, т = 0,1, 2,...,М,
У — ^max
dTm,e,i(r,z,t)
(5)
(6)
z—0
dz
= 0, m = 0,1,2,.
Z — Zm
Здесь и далее Ктах — радиус удаления от траектории иона, а Zmax — глубина, превышающая длину проективного пробега иона. Начальные и граничные условия означают, что в начальный момент времени электронная и ионная подсистемы имеют температуру, равную То, а граница г = Ятах, % = 0, г = Zmax теплоизолирована.
Моделирование динамики фазовых переходов типа «плавление - затвердевание» осуществляется на основе задачи Стефана. Граница раздела фаз определяется условием, что температура вдоль этой границы равна температуре Тшеь
г
фазового перехода, т.е.
Ti(t, £(t, Г, ф, z)) = Tmelt = const > 0, t> 0. (7)
Это соотношение является уравнением для определения £ — положения границы фазового перехода в момент .
Среди методов решения задачи Стефана отметим методы адаптивной сетки, выпрямления фронта и ловли фронта в узел сетки [10]. Однако, упомянутые подходы становятся довольно сложными при переходе к многомерным задачам. Основная трудность при решении этой задачи возникает из-за присутствия движущегося фазового фронта, который в трёхмерном случае представляет собой некоторую поверхность.
В связи с этим для моделирования фазовых переходов используются методы, в которых счёт ведётся без специального выделения границы раздела фаз (разностные схемы сквозного счета). Основную роль в схеме сквозного счета играет принцип сглаживания (размазывания) коэффициентов теплоёмкости Ci(Ti) и теплопроводности Xi(Ti) исходного уравнения. Заметим, что сглаживание коэффициентов производится по Ti и поэтому не зависит, вообще говоря, от того, является задача одномерной или многомерной.
Физическое требование, из которого вытекает граничное условие (7), состоит в том, что при температуре фазового перехода T = Tmeit энергия Н(Ti) как функция температуры испытывает скачок на величину pL, которая называется теплотой (или энтальпией) фазового перехода (здесь р-плотность вещества, L-скрытая теплота).
Следуя [10], введём удельное теплосодержание
Ti
Н(Тг) = J (Ci(x) + pLS(x - Tmelt))dx,
To
C (T ) = / Ci'l(Ti) Ci(±i) I Ci,2(Ti)
X -(T-) = i Xi,l(Ti)
Xi(Ti) = I XiATi)
где ( x) — дельта-функция.
Условия фазового перехода будут учтены, если уравнение (2) заменить уравнением
9Н (Ti) (id ( dTi \ i d2T d2Ti \
~dr = Xi{rd-r{rdi) + ^ a? + ) + g(T - Ti) + Mr,z,^ (8)
Производная Н(Ti) по Ti при Ti — Tmelt обращается в бесконечность, а сама функция Н(Ti) при Ti = Tmeit имеет разрыв первого рода со скачком
Н(Tmelt + 0) - Н(Tmelt - 0) = PL. (9)
Поэтому непосредственное применение разностных схем к уравнению (8) не даёт на практике хороших результатов. Чтобы обеспечить корректное использование стандартных разностных схем для уравнения (8), целесообразно функции Н(Ti), Ci(Ti) и Xi(Ti) подвергнуть сглаживанию.
при То <Ti < Tmelt,
при Tmelt < Ti,
при То <Ti < Tmelt,
при Tmelt < Ti,
Сглаживание осуществляется следующим образом [11]: 1) для энтальпии
CilTi при Ti < Tmelt - ATi, Л
CilTi + (Ti - Tmelt + AT-!)2
AT
Н (Ti) =
(10)
при Tmelt - AT- < Ti < Tmelt--^,
A T A T
C Ti + Е- при Tmelt--^ < T < Tmelt + "y,
- A T 2
C *Ti + El--— {Ti - Tmelt--2 )
при Tmelt + AT/2 < Ti < Tmelt + AT2, „ Ci2Ti + £2 при Ti > Tmelt + AT2;
2) для теплоёмкости
Cii при Ti < Tmelt - ATl,
A T
Cil + Ai(Ti - Tmelt + AT-) при Tmelt - AT1 < Ti < Tmelt--^,
A T A T
Ci(Ti) = { C* при Tmelt--2 ^ Ti < Tmelt + "y,
A T A T
C* - A2 (Ti - Tmelt--при Tmelt + ^ Ti < Tmelt + AT2,
^ Ci2 при Ti > Tmelt + AT2.
В формулах (10) и (11) AT = aATmelt, a < 1, константы e1 и e2 определяются из условия непрерывности функции теплосодержания Н(Ti):
si = Ai(ATi - AT/2)2/2 - (C* - Cii)(Tmelt - AT/2),
£2 = £i - A2(AT2 - AT/2)2/2 + (C* - Ci2)(Tmelt + AT/2),
а константы C*, Ai, —2 определяются по следующим формулам: 1) при Ci2 > Cii
(11)
Н*
2LmeltP - (Ci2 - Cii)(AT2 + AT/2) /(ATmelt + AT); C* =C2i + Н
Ai =
C* C
i
A Ti - A T/2
A2 =
Н*
A T2 - A T/2
2) при Ci2 < Cii
Н*
2LmeltP - (Cii - Ci2)(ATi + AT/2) /(ATmelt + AT); C* = Cii + Н*;
Ai =
Н*
A Ti - A T/2
A2 =
C* C
2
A T2 - A T/2
При численном решении системы уравнений (3)—(4) целесообразно ввести безразмерные переменные, а именно: Те = Те/Т0, Тг = Тг/Т0, г = г/Аг, г = г/Аг, г = Ь/АЪ, где Аг, А г и А1 - единицы измерения расстояния и времени, которые выбраны нами следующим образом: АЬ = 10-13 с и Аг = А г = 10-7м.
*
3. Вычислительная схема
Для численного решения системы уравнений (3)—(6) нами были реализованы три конечно-разностные схемы: стандартная явная схема, схема переменных направлений [12], а также явно-неявная схема, описанная в [9] как «аппроксимация метода матричной факторизации» и имеющая применительно к нашей задаче нижеследующий вид.
Введём равномерную сетку по переменным z, г и t в уравнениях (3)—(4), т.е. положим: ri = i hr (i = 0,1,..., N), Zj = jhz (j = 0,1,... ,L), tk = kht (k = 0,1,... ,K) где hz, hr и ht, соответственно, шаги по переменным z, г и t. Во избежание путаницы индексов температуру электронной подсистемы и температуру решётки в узлах (ri, Zj, tk) обозначим символами Tkj и Tkj (для удобства индекс т опущен). Строим явно-неявную схему следующим образом:
rpk+l _ rpk Г*к i,3 i,3 _ \ °i,3-;- _ Л
h
-1 rpk+l _ rpk+l
1 1i+l,j 1i-l,j , ri+1J
rpk+l _cyrpk+l | rpk+1
+ ri+l,3 21i,3 +ri-l,3 +
+
k k k 1i,j + l 21i,j + 1i,j -1
h?
m
lk+1
2 i, 3
h 2
- atШ -T?.) + Ak+1, (12)
1pk+l _ rpk fjk i,3_ÎJ_ _ д
,
h
h
-, rnk+1 _ rnk+1 rnk+1 _ orpk+l i rnk+1
1 1i+l,j 1i-lJ_ + Ti+l,3 21 i,3 + 1 i-1,3
+
^д k _ r\rpk i ^д k
Ti,j+l 21i,j + Ti,j-l
hl
2
m rpk + 1 hli2 i,j
h г
+
+ -Tkj) + Ak+l. (13)
Начальные условия:
= T0J = 0;i = 0,1, 2,...,N ;j = 0,1, 2,...,L.
Граничные условия при т = 0:
k k k k k k 4 11,j - 1,j - 310,j = 0; 41N-1,j - 1N-2,j - 31N,j = 0;
л^рк rpk nrnk _ ri. лфk rpk nrnk _ ri.
4 11,j - 1 ,j - ö10,j = 0; 41N-1,j - 1N-2,j - 31 N,j = 0;
k = 0,1, 2,..., K; j = 0,1, 2,..., L 41k1 - Тгк2 - 3Tk0 = 0; 4ТгкЬ-1 - ТгкЬ-2 - 3T?tL = 0; 41k1 - Тк - 3Tk = 0; 41* 1 - Тк 2 - ЗТ* = 0;
к = 0,1, 2,..., К ; i _ 0,1, 2, и при m = 0,1, 2,..., M :
,N
k k k k r0,j _ rN,j _ 1i,0 _ 1 i,L,
д k д k д k д k 10,j _ 1N,j _ 1i,0 _ 1i,L,
К.
(14)
(15)
г = 0,1,2,..., N; j = 0,1, 2,...,L; k = 0,1, 2,..
Разностные уравнения (12)—(13) являются неявными только по переменной г и остаются явными по переменной z. Поэтому систему (12)—(13) можно решать применением одномерных прогонок по г. Шаг по времени ht выбирается из условия устойчивости:
h < min 2 Сe,i h2
^ 4 Xet i + h2zg.
Г
Сравнительные расчёты показали, что все три упомянутые конечно-разностные схемы обеспечивают согласующиеся результаты. В частности, получены следующие результаты (Нг = 5 х 10-3, Нг = 4 х 10-2, Н4 = 3 х 10-6):
max
max
Te,yav.s (0, 0, t) — Tt
e,yav.ne
.s(0, 0, t)
Te,yav.s(0, 0, t)
3.4 X 10
-4
Te,yav.s (0, 0, t) -Tt
e,s.per.nap.
.(0, 0, t)
Te,yav.s(0, 0, t)
1.8 x 10
-3
max
Te
e, yav.ne
.s(0, 0, t ) -Te
e,s.per.nap
.(0,0, t)
Te,yav.ne.s(0, 0 t)
1Ax 10
3
Здесь Te.
yav.s
численные результаты явной схемы,
Te
Te
'-е,з.рег.пар. численные результаты схемы переменных направлений, 1 е,уау.пе.з — численные результаты схемы (12)—(13) в точке г = 0, г = 0. Видно, что наибольшая разница получается при сравнении схемы переменных направлений с явной схемой, в то время как рассматриваемая здесь схема (12)—(13) занимает «промежуточное» положение. Дальнейшие расчёты выполнены на основе представленной здесь схемы (12)—(13). Этот выбор обусловлен экономичностью схемы и простотой компьютерной, в том числе параллельной, реализации алгоритма. Для описания фазовых переходов в случае отсутствия аксиальной симметрии необходимо решать систему уравнений (3)-(4) с учётом (10)-(11) для всех т = 0,1,2,...,М.
4. Параллельная реализация
Компьютерная программа, реализующая схему (12)—(13), написана на языке СН—Н с использованием технологии MPI для организации параллельных вычислений. Алгоритм распределения вычислений по MPI-процессам следующий.
Пусть Np — количество MPI-процессов, M — число гармоник, L — количество узлов сетки по .
При малом числе параллельных процессов Np ^ M распараллеливание осуществляется только по индексу m: значения m = 0,..., M — 1 распределяются между процессами. Каждый 1-й процесс вычисляет для назначенных ему значений m полные матрицы (Tij), определяющие температуру в электронной и ионной подсистемах. В этом случае отсутствует обмен данными между процессами во время счета. Взаимодействие процессов происходит лишь на стадии сборки и суммирования по m в рамках процедуры сохранения результатов. Оптимальное распределение нагрузки на процессы обеспечивается, если M кратно Np.
Если Np > M, MPI-процессы распределяются сначала между значениями m = 0,..., M — 1 так, что каждому значению m назначается число процессов
NpmK Далее процессам, предоставленным каждому m, назначаются интервалы
Di = (1 L/N(m) + 1, (1 + 1)L/N(m)) с номерами 1(1 = 0,1, 2,..., N(m) — 1). Каждый процесс Р ведёт расчёты в своих интервалах I для значений jmin(P) < j < j max(P). В каждом временном слое соседние процессы обмениваются рассчитанными граничными значениями фрагментов матриц Tij: каждый Р-й процесс передаёт и получает обновлённые значения Tij, соответствующие граничным значениям j интервалов Di, от соседних процессоров Р — 1 и Р+1 (процессы с номерами 0 и NpT^ — 1 передают и получают граничные значения только от процессоров 1
и NpT^ — 2 соответственно). Сборка и суммирование по m осуществляются в рамках процедуры сохранения результатов в конце счета и на некоторых, заранее заданных, промежуточных современных слоях.
Тестирование эффективности параллельной MPI/C+—Н программы проводилось на кластерах ЦИВК (ЛИТ ОИЯИ, Дубна) и К100 (ИПМ РАН, Москва).
Расчёты (табл.1) показывают, что при Мр = 10 по сравнению с расчётом при Мр = 1 происходит ускорение вычислений примерно в 7-8 раз, что в целом соответствует характерным оценкам эффективности распараллеливания сеточных методов решения уравнений в частных производных. Так, например, в [13] при параллельном решении задачи Дирихле для уравнения Пуассона получено ускорение 5 ~ 6.5 (Мр = 8; 16), а в [14] для одномерного уравнения теплопроводности 5 ~ 7 (^р=9). С увеличением Мр > 20 рост ускорения замедляется и в дальнейшем прекращается из-за возрастающего объёма пересылаемых данных при сборке конечных и промежуточных результатов. Поведение 3(Мр) зависит от значений параметров вычислительной схемы. Для приведённых здесь примеров максимальное ускорение составляет 15-17 раз при Мр = 40.
Таблица 1
Время работы (в минутах) МР1-программы, реализующей МТП, при М = 5, N = 500, т = 0.01 х 10-13с в зависимости от числа МР1-процессов Мр и числа узлов Ь по г. Расчёт проведён на кластере К100 (ИПМ РАН, Москва)
L Np = 1 Np = 5 Np = 10 Np = 20 Np = 30 Np = 40
10000 72.6 19.2 10.0 7.2 5.6 5.1
15000 109.3 26.1 16.1 9.3 8.3 7.2
20000 161.3 38.0 21.3 13.0 11.0 9.5
5. Численные результаты: ионный источник
В этом разделе мы демонстрируем эффект учёта фазовых переходов на примере моделирования облучения никелевой мишени ионами урана с энергией 700 МэВ для аксиально симметричного (2Б) и аксиально несимметричного (3Б) случаев. Численное моделирование позволяет получить профили температур электронного газа и кристаллической решётки в любой момент времени, а также динамику изменения этих температур в любой точке образца.
На рис. 1-2 численные результаты представлены для аксиально симметричного случая. Из рис. 1-2 видно, что профили температур в модели с учётом фазовых переходов существенно отличаются от модели, где фазовый переход не учтён. Горизонтальная штриховая линия показывает температуру плавления никеля, которая равна Тше^ = 1725К. На рис. 3 представлена изотермическая поверхность в случае аксиальной симметрии (а) и в отсутствии симметрии (б) в момент времени г = 3.0 х 10-14с
Т|(0,2Д)К
10 15
(а)
(б)
Рис. 1. Временные зависимости температур кристаллической решётки на глубинах г = 0,4, 8 мкм с учётом фазового перехода (а) и без учёта (б)
0 50 100 150 200 250
(а)
50 100 150 200
(б)
Рис. 2. Радиальные профили температур кристаллической решётки на глубинах г = 0,4, 8,10 мкм в момент времени % = 6 х 10-15 с с учётом фазового
перехода (а) и без учёта (б)
(а) (б)
Рис. 3. Изотермическая поверхность в случае аксиальной симметрии (а) и в отсутствии симметрии (б) в момент времени t = 3.0 х 10-14 c (размеры цилиндра: радиус 140А, высота 200000A)
6. 2D моделирование облучения материалов импульсными
пучками ионов
В этом разделе мы рассматриваем тепловые процессы в материалах при облучении импульсными пучками ионов с учётом фазовых переходов и аксиальной симметрии в рамках нелинейной МТП. Как уже отмечено во введении, здесь используется подход, развитый в разделе 1. Поэтому остановимся лишь на отличии этой задачи от предыдущей. Источник выбираем следующим образом:
Ae,i (t ,Г, Z) = frit) f2(r)S-mel,ph(z),
Где fUt) = — = ^maX 1 - eXP(-«ltf jj, = _1_ Здесь
Ze Ze 1 + exp(«2(i - to))1 1 + exp(«3(r - ro))'
j(t) — плотность тока пучка, f2 (г) определяет распределение пучка по радиусу, параметры jmax, Z,a1,a2,a3 задаются из физических соображений. Переход к безразмерным переменным осуществляется следующим образом: Te = Te/To, % = Ti/To, г = r/Ar, z = z/Az, i = t/At, где At = 10-8c и Ar = Az = 10-5 м.
Далее в качестве материала мишени выбираем железо. Нелинейные параметры для этого случая имеют следующий вид [2]:
( 320 + 0.43Ti + 1.4 х 10-5T?; 300К < % < 1073К, d(Ti) = < 790 + 5.4 х 10-3Tf; 1073К < Tmeit, Tmeit = 1809К, ( 800; Tt > Tmeit,
124 - 0.17Ti + 8.8 х 10-5Т? - 1.3 х 10-8т3; 100К < т < Ттеь,
Чтг)
33; Ti > Tmelt.
В качестве численного метода решения выбираем следующую консервативную схему:
Тк+1 _ тк
çk г,з г,з
,
ht
\к тк+1 _ тк+1 i / \к i \к тк+1 _ тк+1 лг,з тг+1,з тг-1,з + 1 лг+1,з + лг,з тг+1,з тг,з
ihр 2hр h р \ 2 h ^
\к i \к тк+1 _ тк+1 \ i /\к |\к т^к _ т^к
Аг,з + лг-1,з тг,з тг-1,з | + _ лг,з+1 + лг,^ тг,з+1 тг,з 2 hP hz \ 2 hz
\к i \к тк _тк
- дк,э(тЪ ) + Ак+1, (16)
грк+1 _ грк
^к г,з г,з
,
h
\к трк+1 _ трк+1 / \к I \к трк+1 _ трк+1
лг,]тг+1,з тг-1,з , 1 лг+1,з + лг,зтг+1,з тг,з
h p
\ к i \ к тпк+1 _ трк+1
лг i + лг-1,з тг,3
,
г-1,з 1,1 / лг,з+1
) hz V
\ к \ к \ к
1 I лг,i+i + лг,] тг,з+1
\ к - тг,з
2
\к i \к тк _тк
лг,з + лг,з-1 тг,з тг,з-1
+9к?стк? -тк,)+Ак,+^ (17)
Начальные условия:
грО _ грО _ -1
тг,з = т г,з = 1,
г = 0,1, 2,..., N, j = 0,1, 2,...,L.
(18)
Граничные условия:
4т1к_т2^. - зток?-= 0, т^ = 1, 4тк<1-тк<1- 3т^ = 0, т
0,
LN,j
1, - т2,
о,
-N,j
к = 0,1, 2,..., К ; j = 0,1,2,... ,L;
4ткл-тк2 - 3тгко = 0, тгкь = 1, 4тгкл-тгк2 - 3тко = 0, ткь = 1
к = 0, 1, 2, . . . , К, = 0, 1, 2, . . . , N.
(19)
2
z
2
p
h
2
z
2
z
к
1
Эта схема, как и схема (12)-(15), является неявной только по переменной г и остаётся явной по переменной х. Поэтому её решение также сводится к применению одномерных прогонок по г. Как и схема (12)-(15), представленная схема является условно устойчивой.
Параллельный алгоритм является частным случаем схемы, изложенной в разделе 3. Распределение вычислений осуществляется только по индексу ].
Пусть, как и прежде, Ир — количество МР1-процессов, Ь — число узлов сетки по Вводим интервалы индексов^: Б = (IЬ/Ир+1, (¿+1)Ь//Ур); I = 0,1, 2,..., Np— 1. Каждый процесс с номером I вычисляет решения системы (16)-(19) для индексов ], принадлежащих интервалу и каждый процессор I передаёт и получает части массивов искомых решений, принадлежащих граничным значениям интервалов от соседних процессов I — 1 и ¿ + 1 на каждом временном слое (процессы с номерами 0 и Np — 1 передают и получают граничные значения только от процессоров 1 и Np — 2 соответственно).
Рис. 4-6 демонстрируют результаты моделирования тепловых процессов в образце железа под действием импульсного пучка ионов углерода (-2=6), радиус г0 = 3мкм с энергией Е0 = 300кэВ, плотностью тока = 1000Л/сш2 и продолжительностью т = 3 х 10-7с.
3000 2500 2000 1500 1000 500 0
Т,(0,гД)К
%
8 12 16 20
(а)
(б)
Рис. 4. Профили температуры кристаллической решётки по направлению г в центре пучка г — п -—-----™-----------+ — 1 ° 1П-7~ + — ° л -
гз = 3.6 х 10-7с
0 при разных временах ^ = 1.2 х 10 с, = 2.4 х 10 'с > х 10-7с, ÍБ = 6 х 10-7с. (а) и температурные профи направлению г на поверхности г = 0 при тех же временах (б)
(а)
(б)
Рис. 5. Временная зависимость температуры кристаллической решётки на
разных расстояниях от центра пучка г 1 =0 мкм, г2 = 2 мкм, г3 = 4, г4 = 6, гб = 8 мкм на поверхности г = 0 (а) и на разных глубинах г1 =0 мкм, х2 = 4
г3 = 8 мкм, = 12 мкм, ¿б = 16 мкм в центре пучка г = 0 (б)
(а)
(б)
Рис. 6. Динамика изменения радиуса проплава на поверхности г = 0 (а) и глубина проплава в центре пучка г = 0 (б)
На рис. 4 представлены профили температуры кристаллической решётки по г в центре пучка г = 0 (рис. 4а) и температурные профили по г на поверхности г = 0 (рис. 4б) при временах Ь1 = 1.2 х 10-7с, ¿2 = 2.4 х 10-7 с, ¿3 = 3.6 х 10-7 с, ¿4 = 4.8 х 10-7с, ¿5 = 6 х 10-7с.
На рис. 5 приведены временные зависимости температуры кристаллической решётки на разных расстояниях от центра пучка Г1 =0 мкм, Г2 =2 мкм, гз = 4 мкм, г4 = 6 мкм, Г5 = 8 мкм на поверхности г = 0 (рис. 5а) и на разных глубинах
=0 мкм, ^2 =4 мкм, ^з = 8 мкм, ^4 = 12 мкм, ^5 = 16 мкм в центре пучка
г = 0 (рис. 5б). Видно, что продолжительность процесса плавления существенно превышает время обратного процесса затвердевания.
На рис. 6 показана динамика изменения радиуса области плавления на плоскости г = 0 и глубина проплава от оси пучка г = 0.
7. Обсуждение полученных результатов и выводы
Наиболее важным и интересным результатом прохождения высокоэнергетических ионов в твёрдых телах является формирование специфических пространственно распространённых макродефектов - треков. Треком тяжёлой заряженной частицы принято называть сильно деструктированную область вокруг траектории тяжёлого иона в материале, созданную за счёт температурных эффектов, вызванных ионизационными потерями энергии, приводящих к расплавлению с последующей возможной частичной рекристаллизацией или аморфизацией этой области. В ряде материалов, чаще всего в изоляторах, таким способом целенаправленно создают фильтры различного назначения (микродиафрагмы, трековые мембраны и т.д. [15]). Заметим, что, несмотря на довольно длительный период исследований, условия образования треков тяжёлых заряженных частиц высоких энергий в твёрдых материалах неясны до сих пор [16]. В рамках МТП параметры треков принято оценивать как размеры областей, где под действием ионного облучения температура мишени превышает температуру плавления материала мишени.
Экспериментальные результаты облучения для некоторых материалов и сплавов приведены в работах [5,17-19]. В частности, при облучении сплава никеля N1X1 ионами урана с энергией 760МэВ (Т0 = 90К, 5;пе1(^) ~ 57кэВ/нм) диаметр трека составляет 12 нм [17] и при облучении №2гз ионами свинца с энергией 700МэВ ( Т0 = 90К, 5;пе1(^) « 48кэВ/нм) средний диаметр трека 8 нм [18]. В [5] приведены также результаты сравнения экспериментальных исследований облучения различных сплавов железа ионами урана с энергией 8.2 МэВ/нуклон (2 ГэВ) с результатами численного исследования в рамках одномерной МТП. Диаметры трека для сплавов железа, облучаемых ионами урана с энергетическими потерями в диапазоне 5^1(2) ~ 60 — 80кэВ/нм, составляют 8-11 нм. В [19] при облучении хромоникелевой стали ионами ксенона с энергией 124 МэВ диаметр трека составил 10 нм.
В нашей работе представлена параллельная реализация алгоритма для численного решения уравнений МТП, модифицированной для учёта фазовых переходов в материалах под действием ионного и импульсного облучения. Представлены методические результаты, демонстрирующие ускорение вычислений на многопроцессорных системах. Проведённый анализ эффекта учёта фазовых переходов в рамках МТП показывает:
1. При учёте фазовых переходов температура мишени в кристаллических решётках существенно ниже по сравнению с моделью без фазового перехода.
2. При облучении мишени тяжёлыми ионами размеры области, где происходит плавление, следующие: диаметр dшax ~ 15нм, глубина гшах ~ 1.337 х 104нм, в то время как аналогичная область для модели без фазового перехода: йшах ~ 23нм, 2шах ~ 1.4 х 104нм. В этих областях могут происходить структурные изменения (треки, аморфизация, дефекты, разломы) в облучаемых материалах, что может привести к изменению их физических свойств. Полученные результаты не противоречат вышеприведённым экспериментальным данным (характерный диаметр трека 8-12 нм), причём согласие с ними улучшается при учёте фазовых переходов.
3. При моделировании облучения мишени импульсными пучками ионов показано, что длительность процесса плавления на порядок больше, чем длительность обратного процесса - процесса затвердевания, после выключения источника.
4. При облучении мишени импульсным пучком размеры области, где происходит плавление, следующие: диаметр д,шах ^6 мкм, глубина ^шах ~1.6 мкм.
Таким образом, можно сделать вывод, что с учётом фазовых переходов полученные нами оценки размеров областей, где могут происходить процессы плав-
ления, лучше согласуются с приведёнными выше экспериментальными оценками характерных размеров треков по сравнению с моделью без учёта фазовых переходов.
Литература
1. Каганов М. И., Лифшиц И. М., Танатаров Л. В. Релаксация между электронами и решёткой // ЖЭТФ. — 1956. — № 2(8). — С. 232-237. [Kaganov M.I., Lifschitz I.M., Tanatarov L. V. Relaxation Between Electrons and the Lattice // JETF. — 1956. — No 2(8). — P. 232-237. ]
2. Wang Z. G., Dufour C, Paumier E. The Se Sensitivity of Metals under Swift-Heavy-Ion Irradiation: a Transient Thermal Process // J. Phys.: Condens. Matter. — 1994. — Vol. 6, No 34. — Pp. 6733-6750.
3. Toulemonde M., Dufour C., Meftah A. Transient Thermal Processes In Heavy Ion Irradiation Of Crystalline Inorganic Insulators // Nuclear Instruments and Methods in Physics Research Section B. — 2000. — Vol. 166-167. — Pp. 903912.
4. Yavlinskii Y. Track Formation In Amorphous Metals Under Swift Heavy Ion Bombardment // Nuclear Instruments and Methods in Physics Research Section B. — 1998. — Vol. 146, No 1-4. — Pp. 142-146.
5. Morphology of Swift Heavy Ion Tracks in Metallic Glasses / M. D. Rodriguez, B. Afra, C. Trautmann et al. // Journal of Non-Crystalline Solids. — 2012. — Pp. 571-576.
6. Распыление твёрдых тел под действием тяжёлых ионов и температурные эффекты в электронной и решёточной подсистемах / И. В. Амирханов, А. Ю. Дидык, И. В. Пузынин и др. // Физика элементарных частиц и атомного ядра. — 2006. — Т. 37, № 6. — С. 1592-1644. [Sputtering of Solids by Heavy Ions and Temperature Effects in the Electronic and Lattice Subsystems / I. V. Amikhanov, A. Yu. Didyk, I. V. Puzynin et al. // Phizika elementarnih chastic i atomnogo yadra. — 2006. — Vol. 37, No 6. — P. 1592-1644. ]
7. Применение модели термического пика для объяснения изменений структуры поверхности высокоориентированного пиролитического графита при облучении быстрыми ионами 86Kr и 209Bi с высокими ионизационными потерями энергии / И. В. Амирханов, А. Ю. Дидык, Д. З. Музафаров и др. // Поверхность. — 2008. — № 5. — С. 3-12. [Application of Thermal Spike Model for Description of Surface Structure Changes at Highly-Oriented Pyrolytic Graphite under Its Irradiation by High-Speed 86Kr and 209Bi Heavy Ions with High Inelastic Energy Loss / I. V. Amikhanov, A. Yu. Didyk, D.Z. Muzafarov et al. // Poverhnost. — 2008. — No 5. — P. 3-12. ]
8. Численное исследование фазовых переходов, возникающих в металлах под действием импульсных пучков ионов в рамках модели термического пика. / И. В. Амирханов, А. Ю. Дидык, И. В. Пузынин и др. // Поверхность. — 2013. — № 5. — С. 1-6. [Numerical Investigation of Phase Transitions Arising in Metals under an Action of Pulsed Ion Beams in the Frame of Thermal Spike Model / I. V. Amikhanov, A.Yu. Didyk, I. V. Puzynin et al. // Poverhnost. — 2013. — No 5. — P. 1-6. ]
9. Яненко Н. Н. Метод дробных шагов решения многомерных задач математической физики. — М.: Наука-Новосибирск, 1967. — 197 с. [Yanenko N.N. The Method of Fractional Steps for Solution of Mathematical Physics Problems // Moscow: Nauka-Novosibirsk. — 1967. — 197 p. ]
10. Самарский А. А., Вабищевич П. Н. Вычислительная теплопередача. — M.: Едиториал УРСС, 2003. [Samarskiy A. A., Vabishevich P. N. Computational Heat Transfer // Moscow: Editirial URSS. — 2003. ]
11. Численное моделирование динамики температурных полей на плоских мишенях при нестационарном интенсивном лазерном воздействии / М. П. Гала-нин, И. С. Ерхов, Е. Ю. Локтионов и др. // Препринт ИПМ им. М.В. Келдыша РАН. — 2008. — № 61. [Numerical Modeling of Temperature Fields on a Flat Target at Unsteady Intense Laser Pulses / M. P. Galanin, I. S. Erhov, E. Yu. Loktionov et al. // Preprint IPM im. M. V. Keldisha RAN. — 2008. — No 61. ]
12. Самарский А. А., Гулин А. В. Численные методы. — М.: Наука, 1989. — 432 с. [Samarskiy A. A., Gulin A. V. Numerical Methods // Moscow: Nauka. — 1989. — 432 p. ]
13. Гергель В. П. Высокопроизводительные вычисления для многопроцессорных многоядерных систем. — М.: Изд-во МГУ, (серия «Суперкомпьютерное образование»), 2010. [Gergel' V.P. High-Performance Computing for Multi-Core Systems // Moscow: Izd-vo MGU. — 2010. ]
14. Практикум по методам параллельных вычислений / А. В. Старченко, Е. А. Данилкин, В. А. Лаева, С. А. Проханов. — М.: Изд-во МГУ, (серия «Суперкомпьютерное образование»), 2010. [Practical Work on Methods for Parallel Computing / A.B. Starchenko, E. A. Danilkin, V. A. Laeva, S.A. Prokhanov // Moscow: Izd-vo MGU. — 2010. ]
15. Использование ускорителей тяжёлых ионов для изготовления ядерных мембран / Г. Н. Флеров, П. Ю. Апель, А. Ю. Дидык и др. // Атомная энергия. — 1989. — Т. 67, № 4. — С. 274-280. [The Use of Heavy-Ion Accelerators for Producing Nuclear Membranes / G.N. Flerov, P. Yu. Apel, A. Yu. Didyk et al. // Atomnaya energiya. — 1989. — Vol. 67, No 4. — P. 274-280. ]
16. Fink D., Chadderton L. T. Ion-Solid Interactions: Current Status, New Perspectives // Rad. Eff. & Def. in Solids. — 2005. — Vol. 160, No 3-4. — Pp. 67-83.
17. Microstructural Modifications Induced by Swift Ions in the NiTi Intermetallic Compound / A. Barbu, A. Dunlop, A. Hardouin et al. // Nucl. Instrum. Meth. — 1998. — No 145.
18. Latent Tracks do Exist Metallic Materials / A. Barbu, A. Dunlop, D. Lesueur, R. S. Averback // Europhys. Lett. — 1991. — No 15. — Pp. 37-42.
19. Дидык А. Ю. Радиационное воздействие тяжёлых ионов на хромоникелевую сталь при высоких температурах // Известия РАН. Металлы. — 1995. — № 3. — С. 128-135. [Didyk A. Yu. Radiation Effect of Heavy Ions on Chromium-Nickel Steel at High Temperatures // Izvestiya RAN. Metalli. — 1995. — No 3. — P. 128135. ]
UDC 519.624.3
MPI Implementation of the 2D and 3D Simulation of Phase Transitions in Materials Irradiated by Heavy Ion Beams within the Thermal Spike Model
I.V. Amirkhanov, E. V. Zemlyanaya, N.R. Sarker, I. S. Sarkhadov, Z.K. Tukhliev, Z. A. Sharipov
Laboratory of Information Technologies Joint Institute for Nuclear Research 6, Joliot-Curie str., Dubna, Moscow region, 141980, Russia
We present the MPI-based implementation of the method of 2D and 3D calculations of the evolution of temperature fields and the phase transitions dynamics in materials irradiated by high-energy heavy ion and by pulsed ion beams. We utilize a modified thermal spike model based on a system of coupled heat conductivity equation describing thermal processes in the electron and ion subsystems of the target sample. Such equations are numerically solved in the cylindrical coordinate system in axially symmetric (2D) and axially nonsymmetrical (3D) cases. The dynamics of phase transitions is realized on the basis of Stephan's problem in the framework of the enthalpy approach. The mathematical formulation of the problem is given;
a numerical scheme is described; a parallel algorithm is presented. The numerical results confirm the efficiency of our approach and the corresponding MPI/C++ computer code. It is shown that the results of numerical simulations are in agreement with experimental estimations of the track sizes which appear in the target samples exposed to heavy ion.
Key words and phrases: modeling, numerical methods, phase transitions, thermal spike model, parallel algorithm.
