CC BY
49
2010 Дискретные модели реальных процессов №3(9)
ДИСКРЕТНЫЕ МОДЕЛИ РЕАЛЬНЫХ ПРОЦЕССОВ
УДК 621.391.1:004.7
КУМУЛЯТИВНЫЙ СИНТЕЗ: КЛЕТОЧНО-АВТОМАТНАЯ МОДЕЛЬ ПРОЦЕССА ОБРАЗОВАНИЯ ПОКРЫТИЯ, НАНОСИМОГО НА МИШЕНЬ С ПОМОЩЬЮ КУМУЛЯТИВНОГО ПОТОКА ЧАСТИЦ1
О. Л. Бандман*, С. А. Громилов**, С. А. Кинеловский***
* Институт вычислительной математики и математической геофизики СО РАН,
**Институт неорганической химии им. А. В. Николаева СО РАН,
***Институт гидродинамики им. М. А. Лаврентьева СО РАН, г. Новосибирск, Россия
E-mail: bandman@ssd.sscc.ru, grom@niic.ru, skin@hydro.nsc.ru
Предложена дискретная математическая модель стохастического типа, предназначенная для компьютерного моделирования процесса образования покрытия на мишени при обработке её высокоскоростной порошковой струей. Модель относится к классу асинхронных вероятностных клеточных автоматов. Параметры модели определены на основе результатов натурных экспериментов, а также исходя из известных сведений для аналогичных процессов. Результаты моделирования процесса образования карбида на вольфрамовой мишени при обработке её углеродистой порошковой струей показали, что применение модели может повысить эффективность исследований, которые проводятся в рамках создания методов кумулятивного синтеза новых материалов и структур в Институте гидродинамики СО РАН.
Ключевые слова: компьютерное моделирование, асинхронный вероятностный клеточный автомат, кумулятивная струя, твердое покрытие, карбид вольфрама.
Введение
К настоящему времени существующие методы и подходы к синтезу новых соединений, основанные на ударно-волновом нагружении и взрывном компактировании порошковых смесей, в значительной степени исчерпали свои возможности. Это связано прежде всего с ограничениями по диапазонам давлений, температур и массовых скоростей взаимодействующих частиц. Эти диапазоны удается существенно расширить за счет использования специальных кумулятивных зарядов. При изготовлении облицовки кумулятивного заряда из высокопористого материала, в частности из порошковых смесей, реализуемые при обжатии пористых облицовок более высокие, чем при ударно-волновом нагружении, уровни давления и температуры открывают перспективы как для осуществления фазовых переходов в материале, так и для синтеза новых неравновесных структур, отличных от исходного материала [1]. При этом вместо компактной кумулятивной струи образуется разуплотненный поток частиц материала
хРабота поддержана Программой фундаментальных исследований Президиума РАН №2-6 (2010) и Сибирским отделением РАН, Интеграционный проект 32 (2010).
облицовки, который может быть использован для нанесения покрытия на преграду-подложку (мишень). Кроме того, при взаимодействии кумулятивного потока частиц с поверхностью мишени могут образовываться новые фазы и соединения [1]. Совокупность перечисленных результатов дает основание называть данное технологическое направление кумулятивным синтезом. Из изложенного выше следует, что процесс кумулятивного синтеза имеет две стадии, на которых осуществляется собственно синтез:
1) в высокоскоростном потоке частиц (струе) и 2) на поверхности мишени. Данная работа относится к исследованию процессов, происходящих на второй стадии, а именно при образовании покрытия на мишени, обрабатываемой высокоскоростной порошковой струей.
Теоретическая часть исследований включает в себя анализ физико-химических процессов на базе уравнений состояния вещества, а также кинетики фазовых превращений на основе методов молекулярной динамики [2]. Экспериментальная часть состоит из серии специальных экспериментов по взаимодействию высокоскоростного потока частиц с металлической преградой и последующего проведения рентгенофазового анализа образцов полученных покрытий. Результаты таких экспериментов позволяют оценить зависимости структуры покрытия от различных физических параметров процесса. В проведённых экспериментах было осуществлено взаимодействие потока частиц, полученного из 30-градусной конической облицовки из порошка графита насыпной плотности, с металлической вольфрамовой пластиной.
Поскольку организация и проведение экспериментов требуют больших затрат времени и средств, повышение эффективности исследований путем имитационного компьютерного моделирования чрезвычайно важно. Имитационное моделирование позволяет проследить процесс образования покрытия во времени и в пространстве путем наблюдения его в замедленном режиме на мониторе компьютера. Методы компьютерного моделирования, основанные на дифференциальных уравнениях в частных производных, здесь оказываются неприменимыми из-за дискретного характера изменений состояний (мгновенных химических превращений и фазовых переходов). Методы молекулярной динамики полезны для определения количественных оценок процесса механохимического синтеза прекурсоров [2], но не позволяют наблюдать картины процесса на всем пространстве в целом. Наиболее подходящими являются стохастические методы, основанные на непосредственном отображении в компьютере движений и преобразований абстрактных или реальных частиц [3, 4].
Современные компьютеры и суперкомпьютерные системы вполне пригодны для отображения физико-химических взаимодействий на микро- и даже на наноуровнях, при использовании дискретных вероятностных моделей, которые в теории моделирования называются кинетическими асинхронными клеточными автоматами (АКА) [3], а в задачах поверхностной химии — кинетическим методом Монте-Карло [4]. Модель является математическим представлением происходящих в системе событий (адсорбция частицы на поверхность, химические взаимодействия между частицами, фазовые переходы, диффузия по поверхности, проникание внутрь подложки и др.). Каждое событие имитируется изменением состояния покрытия с заданной вероятностью в случайно выбранной точке.
В работе предлагается применить АКА для моделирования процесса образования покрытия, наносимого на металлическую мишень с помощью кумулятивной струи, содержащей порошок из одного или нескольких веществ. В п. 1 дано описание АКА-моде-ли и основанного на ней алгоритма моделирования, а также особенностей программной реализации. В п. 2 приведены результаты проведенного моделирования в усло-
виях, соответствующих указанным выше экспериментам по обработке вольфрамовой мишени высокоскоростной струей из углеродистого порошка. В заключении сделаны выводы о возможностях, целесообразности и перспективах развития методов АКА-моделирования.
1. Описание модели
Поскольку АКА оперирует малыми частицами (размеры — микроны), то, чтобы промоделировать процесс, происходящий на площади в 1мм2, необходим АКА, состоящий из 108 — 109 клеток. Даже при очень простых операциях АКА (несколько обращений к генератору псевдослучайных чисел и несколько булевых операций) обработка такого массива данных требует параллельной реализации на суперкомпьютере, состоящем из нескольких десятков или сотен современных процессоров. Это в настоящее время вполне возможно как с точки зрения вычислительной мощности, так и с точки зрения существования эффективного метода распараллеливания АКА [5]. В случае исследования процессов, в которых носителями энергии являются гранулы порошка (далее называемые частицами) со средним диаметром ~ 10 мкм, для моделирования процессов на участке мишени размерами 1 см2 необходим АКА размерами 1000 х 1000 клеток, что вполне реализуемо на современном персональном компьютере, а для получения картины реаспределения вещества по всей мишени размерами 8 х 8 см2 достаточно иметь 4-ядерный компьютер (например, с процессором типа Intel Core i7).
Задача моделирования процесса образования покрытия на мишени при её бомбардировке частицами порошка состоит в следующем. Мишень представляет собой металлическую пластину, на которую направлена кумулятивная струя, состоящая из порошковой смеси одного или нескольких веществ. При подлете к мишени частицы струи сталкиваются с ней, в результате чего происходят следующие события:
1) адсорбция вещества на поверхность металла;
2) химическая реакция между веществом поверхности и веществом адсорбированной на неё частицы;
3) проникание частицы в мишень на некоторую глубину;
4) перемещение вещества по поверхности покрытия на место с более низким уровнем энергии (диффузия);
5) образование нехимического соединения между веществом поверхности и адсорбированной на него частицы порошка;
6) десорбция (сублимация) вещества с поверхности.
Вероятностный клеточный автомат должен имитировать этот процесс таким образом, чтобы в результате компьютерной реализации модели можно было бы получить полную картину процесса, т. е.
• количество вещества каждого типа в покрытии;
• распределение частиц каждого типа по площади мишени и по слоям покрытия.
Для математического описания алгоритма функционирования АКА удобнее всего
использовать формализм «Алгоритма параллельных подстановок» [6], который, в отличие от классической теории КА, позволяет явно представлять клеточный асинхро-низм (стохастичность изменений как в пространстве моделирования, так и во времени), а также допускает групповые функции изменения состояний, присущие физикохимическим процессам.
Формально АКА определяется тремя понятиями: алфавитом состояний A = {a0, a\,...,an}, множеством имен клеток C = {с} и множеством локальных операторов
0 = {#і(е),... , в<і(о)}. Клетка интерпретируется как пара символов (ак, о), где ак Є А обозначает тип частицы, а о = (і, і, к) —вектор координат дискретного пространства, в котором (і, і) Є М х М — координаты клетки в плоскости мишени, к = 0,... , Н — номер слоя нанесенного покрытия. Номер верхнего (поверхностного) слоя обозначается через к, а Н(і,і) — высота нанесенного слоя в точке (і, і) на плоскости мишени. Таким образом, размер массива равен |С| = Н х М2.
Изменение состояний любой клетки зависит от состояний клеток в её локальном окружении, которое задается шаблоном соседства
Т(о) = {ро(о),рі(о),...^(о)} о = (i,і,k), і,і = 0,...,М,к = 0,...,к(і,і). (1)
где Е {-1,0,11-
Множество локальных операторов 0 = {#1(0),... ,вд (с),... , #п(с)} соответствует множеству событий в моделируемом процессе. Локальный оператор вд (с) Е 0 меняет состояния в некоторых клетках из своего соседства {(уд1, рд1),..., (удг, рдг(с))}, рд1 Е Т(с), дг Е {0,... ,г}, г ^ q, на новые состояния у'д1, равные значениям функций перехода от состояний всех клеток сосдства Т(с), т. е.
где рд — вероятность выполнения вд(с), ь'д1 = /дг (у0,у1, ... ,уд). Так, например, оператор адсорбции меняет состояние поверхностной клетки с = (г,],к(г,])) с V = 0 на V = С. Оператор реакции меняет состояние клетки на результат реакции, если в двух смежных клетках оказываются реактанты, при этом одна из клеток опустошается. Оператор диффузии перемещает частицу в ближайшую клетку, если это перемещение уменьшает энергию системы. Выполнение каждого оператора ограничено условиями, которые вычисляются на основе следующих физических соображений и допущений.
1) Все гранулы порошка имеют одинаковый диаметр В каждой клетке АКА может находиться не более одной гранулы. Физический размер клетки в пространственной решетке равен Пр х Пр. В клеточном автомате все клетки имеют размеры
1 х 1 х 1 при плотном прилегании клеток друг к другу.
2) В струе, падающей на мишень, функции распределения скорости п(г,], ¿) и плотности ¿(г,],Ь) частиц по пространству и времени выражаются либо в виде таблиц, полученных путем компьютерного моделирования струи, либо в виде аппроксимирующих их функций, которые могут быть представлены следующим образом:
При этом принимается р0(с) = с. В кинетических АКА функции рг (с) имеют следующий вид:
Рг(о) = о + (а,в,7), I = 0 ...,д,
г (г, і) = гтахіа ехр(-(Ьг2 + сі)),
(3)
где г (г, ¿) = п(г, ¿) или г (г, ¿) = ¿(г, ¿); г = л](М — г)2 + (М — ] )2 — расстояние от центра мишени; а,Ь,с — константы, различные для п(т, ¿) и ¿(г, ¿).
3) Кинетическая энергия удара частицы о мишень равна
(4)
Эта энергия расходуется на химическую реакцию между веществом упавшей частицы и соприкасающейся с ней на мишени, на деформацию поверхности при проникании
упавшей частицы в глубь мишени, на выделение тепла и, возможно, на образование некоторого нехимического соединения.
4) Вероятность адсорбции частицы на поверхность равна плотности частиц в струе в момент её касания с мишенью
Pads (i,j,t) = d(i,j,t). (5)
5) Химические реакции происходят между частицами, если реактатны находятся в смежных клетках и выполнены условия
Ek(i,j,t) > Ereac(i, j) t') , Treac-low < T(i,j,t) < Treac-high; (6)
где Ereac — энергия активации реакции; Treac-low, Treac-high — границы температурного диапазона, при котором реакция возможна.
6) Перемещение частицы из одной клетки в другую (диффузия) происходит, если после этого перемещения энергия её связей с соседними частицами увеличивается. Поскольку энергию связей между частицами подсчитать практически невозможно, предлагается применить упрощенный вариант диффузии, в котором выполняется только «сглаживание» получаемой поверхности. Иными словами, перемещение частицы происходит с более высокого уровня относительно поверхности покрытия на более низкий, что качественно соответствует принципу уменьшения свободной энергии системы. При этом допустимая разница в уровнях соседних клеток должна быть выбрана исходя из свойств получаемого покрытия. Энергия при диффузии не расходуется.
7) Поскольку химические реакции эндотермические, то можно считать, что когда условия (6) для химической реакции не выполняются или вблизи падающей частицы нет соответствующего реактанта, тогда часть кинетической энергии расходуется на деформацию и нагрев материала:
Ek(i; j,t) Ereac(ii j,t) Edef(i; j,t) + Eheat(ii j,t). (7)
Если не учитывать отвода тепла в течение процесса, то можно считать, что нагревается один слой мишени (толщиной Dp) на участке непосредственно под упавшей на него частицей. При таком допущении повышение температуры можно подсчитать исходя из теплоёмкости материала мишени и массы нагреваемого материала.
8) Энергия деформации поверхности расходуется на проникание в глубь мишени и на пластическую деформацию вблизи падения частицы, т. е.
Edef(i,j,t) Esurf ' Sbr(i,j,t) + Educ ' Vduc(i,j,t), (8)
где Esurf — поверхностная энергия изменяющего свою форму материала; Sbr — изменение площади поверхности; Educ — энергия пластической деформации; Vduc — объём пластически деформируемого материала. Величины энергий Esurf и Educ характеризуют хрупкость и пластичность материала.
2. Реализация программы и результаты моделирования
Модель испытывалась на примере вольфрамовой мишени, обрабатываемой углесодержащей кумулятивной струей. Исходные данные для моделирования соответствуют экспериментам, которые были проведены в Институте неорганической химии и Институте гидродинамики СО РАН: размер мишени 80 х 80 мм2; максимальные скорость и плотность порошка, падающего на мишень, составляют umax = 2000 мс-1, dmax = 0,4;
средний размер гранул d =10 мкм; струя содержит ~ 2,8-109 гранул порошка углерода с массой гранулы т = 1,8 • 10-9 г.
Программная реализация отлаживалась на фрагменте размерами 6 х 6 мм2. При линейном размере модельной клетки I = 10 мкм размер фрагмента составляет 600 х 600 клеток в плоскости мишени при максимальной высоте наносимого покрытия Н(г,]) = 30. Такие размеры АКА позволяют провести экспериментальную отладку программы на персональном компьютере типа РеПлиш IV при выводе на экран монитора моделируемой поверхности на каждой итерации. Более того, если предположить, что выбранный фрагмент находится в центре мишени, то значения скорости и плотности в струе можно считать постоянными и равными своим максимальным значениям. Именно эти значения можно сравнивать с данными дифрактограмм результатов натурных испытаний. При этом предположении и с учетом известных экспериментальных данных определены следующие параметры модели.
1) Кинетическая энергия, приносимая падающей частицей: Е^ = 3,6 • 10-6 Дж.
2) Длительность существования струи Ь = 25 мкс при расстоянии от облицовки до мишени Ь = 50 мм.
3) Вероятность адсорбции Раа8 = ртах = 0,4 равна плотности струи при соприкосновении с мишенью.
4) Вероятности реакций могут быть рассчитаны исходя из соотношений значений энергий активации (энтальпий) реакций карбидов Еа(ШС) и Еа(Ш2С) при заданных условиях (температуре, давлении, количестве избыточного углерода). Однако поскольку, кроме предположительной температуры в центральной части струи Ттах = = 1600-2000 ° С, других данных для локальных условий моделируемого процесса нет, то приходится ориентироваться на литературные данные [7, 8] и данные натурных испытаний авторов, из которых можно заключить, что при заданной температуре и избытке углерода соотношение между количествами получаемых карбидов и остающимся вольфрамом равно ф(Ш)/ф('ШС)/ф('Ш2С) = 2/3/4, что позволяет принять следующие величины вероятностей образования карбидов: Р(Ш2С) = 0,44, Р(ШС) = 0,33.
5) Диффузия выполняется с вероятностью Ращ =1 в тех клетках массива, где разница высоты покрытия в соседних клетках АН ^ 2. Акт диффузии состоит в том, что частица с более высокого уровня перемещается на более низкий, причем если существует несколько возможных вариантов перемещений, то один их них выбирается равновероятно.
6) Из экспериментов известно, что проникание углерода в глубь вольфрамовой ми-
шени происходит не более чем на несколько десятков микрометров. Таким образом, энергия деформации Е&^ рассчитывается по (8) при условии проникания поверхностной клетки на два слоя вниз и вытеснения аморфного вольфрама на поверхность. При ЕзигКШ) = 1,7 Дж-м-2, Еаис(ШС) = 7МПа- Н/м3/2 и £Ьг = 0,9 • 10-9 м2, = 10-12 м3,
рассчитанных по данным из [8] и Т = 1600 ° С, энергия составит Еа^ = 2,2 • 10-9 Дж. Остальная энергия расходуется на повышение температуры в месте соударения. Зная теплоемкость вольфрама cw и принимая нагреваемую массу мишени тщ равной массе вещества в одной клетке, получаем повышение температуры АТ = ЕЬетт/(сщ • тщ). Эти расчетные величины могут быть приняты в качестве ориентировочных, так как они зависят от свойств и получаемых карбидов, и мишени, которые значительно различаются от образца к образцу.
Функционирование АКА происходит в соответствии со следующим стохастическим итерационным алгоритмом.
В начальном состоянии (t = 0) все клетки массива П(0) пусты:
y(i,j,k) е С (v(i,j, 0) = 0,h(i,j) = 0).
Каждая t-я итерация состоит из трех частей:
1. Для всех (i,j) е С, выбираемых в случайном порядке, выполняется следующее:
• с вероятностью Pads = pmax = 0,4 применяется оператор адсорбции;
• если условия (6) для реакций выполнены, то в соответствии с вероятностями
P(W2C) и P(WC) применяется одна из химических реакций: W+C ^ WC или W+W+C ^ W2C;
• если условия для химических реакций не выполнены, то применяются операторы
деформации и нагрева.
2. Для всех (i,j) е С, выбираемых в случайном порядке, выполняется операция диффузии.
3. Через несколько итераций покрытие оказывается сформированным, так как поверхность оказывается покрытой слоем углерода и химических реакций не происходит. Процесс заканчивается полностью, когда иссякает поток порошка.
Программная реализация предусматривает получение следующей информации о процессе.
1) Вывод в файл зависимости процентных соотношений количеств W, WC и W2C в поверхностном слое и во всем покрытии от времени (на каждой итерации).
2) Вывод на экран монитора на каждой итерации состояния поверхностного слоя и разреза мишени, причем каждое состояние клетки (вещество) представлено своим цветом.
На рис. 1 приведены зависимости количеств полученных веществ в поверхностном слое от времени моделирования. Время выражено в количестве итераций, каждая итерация соответствует реальному промежутку времени At = 10-8 с.
На рис. 2 показано распределение веществ по уровням (слоям) покрытия, отсчитываемым от h(i,j) = 0, причем номер слоя указывает на расстояние от нулевого уровня, т. е. уровня поверхности вольфрама до начала процесса. Из рис. 2 видно, что наибольший процент карбидов находится в слоях ниже нулевого, что объясняется тем, что получающийся в результате реакции карбид занимает место участвовавшего в реакции вольфрама.
Моделирование процесса на мишени размерами 80 х 80 мм2 проводилось на компьютере с 4-ядерным процессором типа Intel® Core™ i7. Клеточное пространство размером 8000 х 8000 было разделено на четыре домена размерами 4000 х 4000 каждый. Распараллеливание асинхронного алгоритма клеточно-автоматной эволюции проводилось с применением блочно-синхронного преобразования АКА [7]. При этом учитывалось распределение скорости и плотности частиц в струе по пространству, которые рассчитывались в соответствии с формулой (3). Константы в (3) a = 1,86, с = 0,64, bu = 30 и bd = 12 были вычислены путем подстановки в (3) следующих граничных значений: u(r = 0,t = 2) = umax, d(r = 0,t = 2) = dmax, u(r = 0,t = 20) = 0, d(r = 0,t = 20) = 0, u(r = M/2) = d(r = M/2) = 0 для всех t = 0,..., 20.
Итерации
\Л/ -*-\Л/С —\Л/2С
Рис. 1. Зависимость количества веществ в поверхностном слое от времени
Номер слоя -^-УУС УУ2С
Рис. 2. Распределение веществ в покрытии по слоям
При моделировании применялся тот же алгоритм, что и для центрального фрагмента. Но для каждой точки с координатами (і, j) вычислялись значения и (і, ^') и сІ(і^)
и, в зависимости от них, вероятности
РаОв&Л = ¿(і,і), Рпс(і,Л = 0,4ехр(1,4 - г;2(г,і)), Р№2С = 0,6ехр(-2,4 - г;2(г,і)),
где г’(ї, і) = и (г, і) • Ю_3. Полученные в результате моделирования распределения скорости и плотности струи вдоль её оси показаны на рис. 3.
X
го
Е
х
го
£
Итерации
Члах<«
^тах(^)'0:00^
Рис. 3. Зависимость скорости и плотности в центре струи от времени ее касания с мишенью
На рис. 4 показаны полученные в результате моделирования соотношения плотности веществ в зависимости от расстояния от центра мишени. Время вычисления на компьютере Intel® Core™ i7 составляет 0,5 ч, что убеждает в возможности проводить в приемлемое время моделирование процессов с порошками размером 1 мкм и менее.
120
Л?" лГ ^ &
^ <Г
Итерации
Рис. 4. Распределение веществ по пространству мишени после окончания процесса
Заключение
Представлена дискретная математическая модель стохастического типа, относящаяся к классу асинхронных клеточных автоматов, предназначенная для моделирования процесса образования покрытия на мишени при обработке её высокоскоростной порошковой струей. Основные параметры модели, а именно вероятности адсорбции и химических реакций, скорость, плотность и температура порошковой струп приняты в модели исходя из натурных экспериментов и полученных на их основе дифракто-грамм. Однако часть необходимых данных (поверхностная энергия, энергия пластической деформации), а также допущений (упрощенная диффузия, функция распределения скорости и плотности порошка в струе) были выбраны исходя из известных сведений для аналогичных процессов [7, 8]. На самом деле, для каждого набора порошков и материала мишени они должны быть уточнены путем проведения серии натурных и соответствующих вычислительных экспериментов. Программные реализации,
разработанные для малого фрагмента мишени и для мишени в целом, показали, что она может стать полезным инструментом в исследованиях процессов кумулятивного синтеза.
ЛИТЕРАТУРА
1. Громилов С. А., Кинеловский С. А., Попов Ю. Н., Тришин Ю. А. О возможности физикохимических превращений веществ при кумулятивном нанесении покрытий // Физика горения и взрыва. 1997. Т. 33. №6. С. 127-130.
2. Псахье С. Г., Коростелев С. Ю., Смолин А. Ю. и др. Метод подвижных клеточных автоматов как инструмент физической мезомеханики материалов // Физическая мезоме-ханика. 1998. T. 1. №1. C. 1-15.
3. Бандман О. Л. Клеточно-автоматные модели пространственной динамики // Системная информатика. Новосибирск: СО РАН, 2006. Вып. 10. С. 59-113.
4. Jansen A. P. J. An Introduction to Monte Carlo Simulation Of Surface Reaction // arXiv: cond-mat/0303028v1[cond0-mat.stat-mech]
5. Bandman O. Parallel Simulation of Asynchronous CellularAutomata Evolution // LNCS. 2006. V. 4173. P. 41-48.
6. Achasova S., Bandman O., Markova V., Piskunov S. Parallel Substitution Algorithm: Theory and Application. Singapore: World Scientific, 1994. 180 p.
7. Itaka I., Aoki Y. Quantitative separation of WC from W2C and tungsten, and the conditions of formation of two carbides // Bulletin of Chemical Society. 1982. No. 4. P. 108-114.
8. http://www.wolfram_at/wDeutsch/produkte/carbid/WC_eigenschaft. Common Properties — Tungsten Carbide.
