Создание новых материалов с заданным и управляемым комплексом свойств является приоритетным научным направлением. Для металлических материалов эти свойства во многом определяются характеристиками жидких фаз металла и шлака и закономерностями обменных взаимодействий на границе их раздела. Основой большинства металлургических шлаков, являются многокомпонентные оксидные расплавы, являющиеся важными объектами исследования в физической химии.
Эти расплавы относятся к типу неупорядоченных сильновзаимодействующих ионно-ковалентных полимеризующихся систем, большинство свойств которых является структурочувствительными. Синтез новых металлических материалов может быть эффективно осуществлен лишь на основе знания зависимостей физико-химических свойств оксидов и их структурных особенностей. В данной предметной области имеется большой разрыв между возможностями теории и эксперимента. Этот разрыв с успехом заполняет компьютерный эксперимент, часто оказывающийся единственным способом изучения протекающих процессов и явлений. Поэтому построение методов и моделей для компьютерного моделирования, развивающих теоретический фундамент физической химии и расширяющих класс практически важных результатов численного эксперимента, является крайне актуальным.
Целью являлась разработка методов компьютерного моделирования наноструктуры многокомпонентных оксидных расплавов, интегрирующих подходы полимерной теории и теории графов, с использованием результатов молекулярно-динамического моделирования фазовых траекторий системы многих частиц и их применение в информационно-исследовательской системе «MD-Slag-Melt» для расчета практически важных свойств.
Современные молекулярно-динамические (МД)-пакеты (LAMMPS, HyperChem, XMD и др.) [5] позволяют проводить компьютерные эксперименты для систем, содержащих миллионы частиц, используя высокоскоростные параллельные вычисления. Как правило, в результате получают файлы координат и скоростей, размером в несколько гигабайт, по которым путем статистической обработки строят функции радиального распределения, а также некоторые термодинамические параметры(температура, давление, теплоемкость) и коэффициенты диффузии.
Для реальных нужд существенно важны транспортные свойства, которые являются структурочувствительными. Однако моделирование многокомпонентных оксидов значительно затрудняют специфические особенности (большие значения вязкости, характер электро- и теплопроводности, экспоненциальная зависимость коэффициентов переноса от температуры), которые связаны со структурной неоднородностью среднего порядка и медленно развивающимися процессами, происходящими между полианионными комплексами, объединяющими от нескольких до сотен или тысяч частиц с ионно-ковалентными связями.
Для решения задач моделирования структуры и свойств шлаковых расплавов разработан программный комплекс с удаленным доступом - информационно-исследовательская система (ИИС) «MD-SLAG-MELT»[1]. ИИС содержит систему классов математических методов и моделей, обеспечивающих комплексное моделирование расплава, в том числе модели межчастичного взаимодействия[4], молекулярно-динамическую модель, модели структуры, модели физико-химических свойств.
Для теоретического обоснования и предсказания поведения и свойств исследуемых систем на основе данных МД-эксперимента необходимо расширение класса моделей структуры до моделей, описывающих особенности структуры среднего порядка (наностуктуры). Это требует разработки математического формализма, обобщающего уже существующие модели на неоднородный уровень, совокупности специфических вычислительных методов и алгоритмов высокого уровня сложности.
Ниже описывается формализм разработанного авторами метода для описания структуры многокомпонентной ионно-ковалентной оксидной системы по результатам молекулярно-динамического моделирования, интегрирующего подходы полимерной теории и теории графов.
Метод моделирования неоднородной структуры среднего порядка на основе модели обобщенных дескрипторов в неоднородной графовой сети опирается на модель структуры расплава в терминах неоднородной графовой сети, с поэтапным наращиванием сложности отображения от выделения разнородных графовых вершин, описывающих исследуемые в модельном кубе объекты, до формирования компонент связности и количественной оценки структурочувствительных характеристик. Сначала выделим основные этапы метода, а затем более подробно опишем содержание каждого:
1) построение множеств вершин на основе неоднородного дескриптора;
2) выделение звезд характеристических вершин;
3) формирование связных графов «второго уровня»;
4) поиск минимальных замкнутых контуров в связных графах «второго уровня»;
5) формирование модели структурочувствительных свойств.
Построение множеств вершин на основе неоднородного дескриптора. Изначально, на этапе построения множеств графовых вершин, строится неоднородный дескриптор d, содержащий элементы, отображающие описание структуры в терминах полимерной модели и сведений об МД-эксперименте на теоретико-графовую нотацию:
,
где type - тип вершины, kmax - максимально допустимая валентность вершины, kcurr - текущая валентность вершины, i - порядковый номер вершины, r - радиус координационной сферы вершины (или характерное расстояние между вершинами определенного типа), kf - скалирующий множитель, s - агрегатор, описывающий типы смежных вершин.
Множество объектов модельного куба разбивается на разнотипные подмножества графовых вершин по заданному неоднородному дескриптору. На первом этапе объекты разделяются на 3 подмножества вершин:
- характеристические (структурообразующие)
при
;
- модифицирующие
при
;
- однородные
при
.
Однородные вершины в соответствии с каждым значением максимальной валентности подразделяются на изолированные 
, висячие 
и двухвалентные 
.
К подмножеству вершин дополнительно применяется агрегатор s, описывающий тип и валентность смежной вершины. При 
, выделяется множество однородных конечных двухвалентных вершин 
.
Модифицирующие вершины подразделяются на два типа в соответствии с валентностью: при kmax = 1 - висячие модифицирующие вершины 
; при kmax > 1 - изолированные модифицирующие вершины 
.
Этап построения вершин определяет степень неоднородности исследуемой структуры и является основой для её дальнейшего исследования методами теории графов.
Выделение звезд характеристических вершин. Следующим этапом является выделение звезд для всех характеристических 
и модифицирующих вершин 
. Звезды являются структурными единицами среднего уровня и описывают связи между разнотипными множествами вершин графа.
В модельном кубе для каждой характеристической (и модифицирующей) вершины выделяется звезда 
, на концах которой могут находиться только вершины из множества однородных:
,
причем количество смежных вершин не должно превышать максимальную валентность 
, все они должны находиться в пределах радиуса первой координационной сферы r. Звезды являются базовыми компонентами построения связных графов.
Формирование связных графов «второго уровня». Одно из заданных условий метода предполагает, что звезды вершин могут объединяться в более крупные обрзования - связные графы. С этой целью на множестве звезд вводится теоретико-графовая структура - граф «второго уровня» G, вершинами которого являются звезды вершин.
.
Две звезды 
и 
являются связными при совпадении однородных двухвалентных вершин
.
На графе «второго уровня» происходит обобщение множеств разнотипных вершин и поиск компонент связности, характеризующих структурную сеть модельного куба. На основе выделенных связных компонент проводится исследование разрывов структурной сети.
Поиск минимальных замкнутых контуров в связных графах «второго уровня». Замкнутым контуром в связном графе «второго уровня» 
является минимальная последовательность неповторяющихся связных вершин, начальный и конечный элементы которой совпадают
.
Замкнутые контуры описывает плоскостные кольца в графовой сети и являются одним из показателей неоднородности структуры.
Модель структурочувствительных свойств. Комплекс построенных характеристик - распределение звезд вершин, связных графов и замкнутых маршрутов в них позволяет построить систему дескрипторов, отображающих физико-химические свойства расплавов. В модели структурочувствительных свойств выделено 2 блока параметров: функции распределения графовых вершин и скалярные структурные характеристики.
Функции распределения вершин связных графов по элементам дескриптора dG.
,
где 
- доля связных графов с одинаковым числовым значением дескриптора; k - количество конфигураций; 
- количество связных графов с одинаковым числовым значением дескриптора; τi - время жизни i-го связного графа.
Скалярные структурные характеристики
1. Доли однородных вершин в модельной системе:
, Nkonf - полное число конфигураций.
2. Степень полимеризации системы 
, константа полимеризации - 
, где X(On) - доли разнотипных однородных вершин.
3. Энергии смешения оксидов, образующих расплав, изменения энтальпии и энтропии системы, активности участвующих оксидов: для сеткообразующего
,
для модифицирующего 
,где 
, 
- соответственно относительные числа характеристических звезд вершин, содержащих k двухвалентных однородных вершин и изолированных однородных вершин к общему количеству графовых вершин в системе.
4. Доля замкнутых контуров в связном графе
,
где l - количество связных графов.
Алгоритмизация метода. Для поиска звезд вершин разработан алгоритм выделения областей Шлегеля, существенно ускоряющий исследование особенностей графовой сети, за счет поиска смежных вершин для каждой структурообразующей вершины только в областях Шлегеля. Каждая вершина описывается номером, типом и координатами: 
. Модельный куб разбит на нумерованные ячейки, ассоциированные с номерами однородных вершин, в них расположенных, с учетом граничных условий:
,
где xj, yj, zj - координаты вершины; 
- номер ячейки.
Поиск идет по области Шлегеля, связанной с ячейкой центральной вершины, а также 26 смежных с ней ячеек.
,
где M6(m) - соседи по граням ячейки, M12(m) - соседи по ребрам ячейки, M8(m) - соседи по вершинам ячейки. Таким образом, можно утверждать, что
- область Шлегеля для вершины 
.
После выделения области Шлегеля для каждой характеристической вершины, осуществляется поиск ближайших соседей, с ограничениями на валентность центральной вершины и радиус координационной сферы. Полученные звезды описываются номерами вершин.
В таблице приведены результаты тестирования алгоритмов поиска звезд методами простого перебора и выделением областей Шлегеля. Сложность алгоритма возрастает с увеличением количества неоднородных дескрипторов и составляет:
,
где k - количество конфигураций; N - количество вершин, N0 - количество однородных вершин; L - длина модельного куба.
Результаты тестирования
|
Количество частиц |
Простой |
Поиск |
|
500 |
0,08 |
0,11 |
|
1 000 |
0,15 |
0,18 |
|
2 000 |
0,58 |
0,38 |
|
4 000 |
2,3 |
0,83 |
|
5 000 |
3,6 |
0,97 |
|
10 000 |
15,7 |
2,4 |
|
20 000 |
89,5 |
6,3 |
|
40 000 |
594,9 |
19,8 |
Следующим этапом является выделение особенностей структуры среднего порядка, что обеспечивается методом исследования связности графов «второго уровня».
Алгоритм выделения кратчайших замкнутых маршрутов. На основе полученных характеристик о наличии связных графов, каждый из которых описывается списками смежности вершин и временем жизни, осуществляется поиск минимальных замкнутых маршрутов, реализацией рекурсивного алгоритма поиска в ширину (BSF) на графе «второго уровня».
Пусть 
- обход при поиске в ширину, а D(gi, j) - уровни обхода графа относительно вершины 
- вершина, принадлежащая D(gi, j) уровню обхода графа. Тогда при обходе графа по j-й вершине условиями существования кольца являются: 
- при условии четного количества элементов в кольце, 
- при условии нечетного количества элементов.
При этом не учитываются кольца, являющиеся соединением двух и более соседних колец, что позволяет значительно сократить время поиска минимальных колец, исключая более крупные образования. Разработанный алгоритм разбиения трехмерного пространства в модельном кубе и поиска ближайших соседей,
Разработанные методы и алгоритмы реализованы в ИИС «MD-Slasg-Melt»[1]. Ниже приведены результаты моделирования системы SiO2-Na2O в диапазоне шести составов(0,4-0,6; 0,43-0,67; 0,5-0,5; 0,63-0,37; 0,74-0,26; 0,8-0,2), которая играет важную роль в металлургической и стекольной промышленности. Температуры моделирования превышают на 30-50 °С температуру плавления.
Метод построения графовых вершин позволяет рассчитать характеристики структуры ближнего порядка, в частности координационные числа (к.ч.) (рис. 1а) и валентные углы (рис. 1б), из которых следует, что к.ч. кремния практически не зависят от состава, а к.ч. натрия увеличиваются при уменьшении мольной доли оксида натрия в составе.
а б
Рис. 1. Расчет характеристик структуры ближнего порядка:
а - координационные числа для Si и Na;
б - распределение валентных углов для O-Si-O и Si-O-Si
Углы O-Si-O в среднем равны 109°, что соответствует sp3-гибридизации. Углы Si-O-Si, т.е. углы между соседними структурными группами, лежат в пределах от 120 до 180°.
С помощью метода моделирования наноструктуры для всех составов получены функции распределения комплексных анионов, а также исследована их протяженная структура, компоненты связности, их время жизни и скалярные параметры.
Данные по распределению разных типов кислорода по составам говорят о том, что концентрация свободного кислорода снижается до нуля при NSiO2 = 0,7. Характер изменения содержания концевого и мостикового кислорода противоположен и свидетельствует об увеличении степени полимеризации системы при увеличении мольной доли оксида кремния, вплоть до непрерывной сетки. Этот процесс отражается на рис. 2, где приведены результаты распределений полианионных комплексов по типам и времени жизни.
а б
Рис. 2. Результаты обработки распределений вершин на всех конфигурациях:
а - доли комплексов; б - время жизни комплексов
Процесс полимеризации выражается в том, что базовые структурные группировки SiO4 объединяются в полианионы нарастающей степени сложности. В разработанной модели структуры учтена особенность атома натрия, который присоединяется к концевым атомам кислорода в полиаонинах.
Образующиеся полианионные комплексы постоянно обновляются, так как время жизни даже крупных группировок (более 50 атомов кремния) в средней области составов не превышает 30 % общего времени моделирования. Это связано с миграцией атомов натрия.
Константа полимеризации уменьшается при увеличении содержания SiO2. При постепенной полимеризации оксида реакция преобразования мостикового и свободного кислорода в два концевых затрудняется, так как миграции ионов внутри сетки практически нет, а количество свободного кислорода уменьшается от состава к составу.
Выводы
В статье рассмотрен формализм методов для компьютерного моделирования структуры среднего порядка(наноструктуры) многокомпонентных оксидных расплавов, интегрирующий подходы полимерной теории и теории графов на основе данных молекулярно-динамического моделирования. Полученные теоретические результаты обладают научной новизной
Кроме того в статье приведены результаты применения описанных методов к моделированию бинарной системы SiO2-Na2O, отражающие особенности процесса полимеризации и обладающие практической значимостью.
Рецензенты:
Ивашко А.Г., д.т.н., профессор, зав. кафедрой информационных систем Института математики и компьютерных наук ГОУ ВПО «Тюменский государственный университет», г. Тюмень;
Копыльцов А.В., д.т.н., профессор, зав. кафедрой информатики ГОУ ВПО РГПУ им. А.И. Герцена, г. Санкт-Петербург.
Работа поступила в редакцию 27.06.2011.