Системы и технологии управления параллельными вычислениями (стр. 9 )

Рис. 8.1. Колебания давления и «средние» при разных температурах в NVT МД.

Рис. 8.2. Отклонения давления и периода от их экспериментальных температурных зависимостей, приведенные к одному масштабу.

Поэтому для решения оптимизационной задачи за наименьшее время мы использовали NVT-моделирование (при постоянстве количества частиц, объёма и температуры) и минимизировали отклонения от эксперименталь­ных значений уравнения состояния по давлению, вместо объема (периода решетки), а NPT-моделирование (при постоянном нулевом давлении вместо постоянства объёма) температурной зависимос­ти периода решетки приме­няли лишь для заключительной оценки качества полученных ПП.

На рис. 8.1 показаны характерные кривые колебаний давления для NVT-моделирования диоксида урана при температурах 300 К (комнатной) и 3100 К (близкой к плавлению) и их усреднение по времени. Видно, что на временах порядка 10 пс «средние» приходят к своим равновесным значениям.

На рис. 8.2 приведены графики отклонения периода и давления от экспериментальных при одинаковом масштабе, видно, что кривая давления (в NVT) воспроизводит все особенности температурной зависимости отклонения периода (в NPT) – знаки функции и первых двух производных. Тем не менее, можно отметить количественные отличия между кривыми, нарастающие с температурой, эти отличия обусловлены нелинейной зависимостью уравнения состояния твердого тела от температуры. Поэтому оценка отклонений периода по давлению на высоких температурах является заниженной (максимальное отличие наблюдается вблизи плавления и достигает ~0.01 Å).

9.3.  Модель и детали реализации

Парные потенциалы искали в форме Борна-Майера [38], при этом для снижения размерности оптимизационной задачи использовали распростра­ненное приближение [39]–[43]:

· у взаимодействия «катион-анион» пренебрегали дисперсионным притяжением, а у «катион-катион» обоими составляющими в сравнении с кулоновским членом;

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

Расчет межчастичных взаимодействий в ПГУ проводили с помощью суммирования Эвальда в приближении «ближайших отражений» [44]:

Здесь r – радиус-вектор, k – вектор обратной решетки, Uij – короткодействующий ПП в форме Борн-Майера с подбираемыми параметрами Aij, Bij, Cij, Q – подбираемый безразмерный множитель зарядов qi и qj, Ke = 14.3998 эВ*A – электростатическая постоянная, а w – параметр метода Эвальда.

Исходная транслируемая ячейка содержала 768 ионов – 4х4х4 элемен­тарные ячейки по 12 частиц, т. к. дисперсионное слагаемое потенциала анион-анион при 3х3х3 (и меньшем числе ячеек) на границе транслируемой области вносит вклад порядка ~10–3 в расчет межчастичных сил и его обрезание на этом расстоянии приводит к постоянному увеличению энергии (разогреву) в системе. Сумма в обратном пространстве бралась по 2196 векторам обратной решетки. При таких параметрах относительная точность вычисления сил, действующих на частицы, составляла: ~10–6.

Динамику системы моделировали методом Рябова [45] с интегратором Верлета [46]:

Здесь rFCC – начальные положения ионов, образующие гранецентриро­ванную кубическую решетку с экспериментальным значением объема V0 (соответствующего экспериментальному периоду решетки); vMB – начальные скорости, соответствующие распределению Максвелла-Больцмана для выбранной температуры T0; M, P, U, K – масса, давление, потенциальная и кинетическая энергии модельного кристалла в момент времени t; ∆t = 5 фемтосекунд – шаг времени интегратора Верлета; g = 2 – безразмерная константа баростата Рябова, определяющая частоту колебаний для радиуса гипертора R (R¢ – его скорость).

Во время релаксации системы к равновесию каждые 0.2 пс (т. е. примерно каждые 2 периода тепловых колебаний ионов кислорода) прово­дилась коррекция скоростей в соответствии с заданной температурой T0 масштабированием на (T/T0)1/2. После релаксации в течение времени усред­нения накапливали средние значения следующих макропараметров [47]:

Здесь Т – температура, kB – постоянная Больцмана, BT – изотермичес­кая сжимаемость, E и H – общая энергия и энтальпия системы,CV и CP – изохорная и изобарная теплоемкости, DS – коэффициент самодиффузии частиц выбранного сорта число которых – S.

Параметры ПП искали в виде решения оптимизационной задачи мини­мизирующего среднеквадратичное отклонение от нуля среднего давления по 4-6 точкам температуры в диапазоне 300~3100К:

Локальную нелинейную многомерную оптимизацию проводили мето­дом Брента [48], который подобно распространенному симплекс-методу Нельдера-Мида [49] не требует производной, но демонстрирует лучшую сходимость. Для глобальной оптимизации начальный набор параметров локального метода определялся случайным отклонением от найденного минимума по разным переменным с учетом их масштаба. В качестве начальных параметров ПП для оптимизационной задачи использовали значения, полученные другими авторами [39]–[43] [50]–[51], а затем – ранее найденные в процессе оптимизации.

В результате, на основе разработанных методики и высокопроизво­дительной (см. табл. 8.2) технологии МД-моделирования на графических процессорах MD-GPU [52]–[55] нам удалось в течение полугода получить наборы парных потенциалов для разных типов оксидного ядерного топлива: UO2, PuO2, MOX.

В данной работе для обсуждения предложенной методики и получен­ных результатов мы решили подробно проанализировать диоксид урана, т. к. для этой системы доступно наибольшее число экспериментальных данных и работ других авторов. Также это соединение в настоящий момент наиболее важно с технологической точки зрения, т. к. на его основе изготовлено более 70% всех использующихся ТВЭЛов. Параметры парных потенциалов, полу­ченные нами для этого соединения, приведены в табл. 8.3.

Таблица 8.2

Затраты времени на решение оптимизационной задачи

Технология расчетов	1 шаг МД (768 ионов)	шагов на 1 точку	4-6 точек Т на 1 итерацию поиска	100000 итера­ций поиска на 1 ПП
MD-GPU	0.005 сек	4-12 сек	16-72 сек	18-83 суток
MD CPU	0.21 сек	168-504 сек	11-50 мин	2.1-10 лет

CPU: AMD Athlon64 3200+. GPU: ATI Radeon X1900XT 512MB.

Таблица 8.3

Параметры полученных эмпирических парных потенциалов
для диоксида урана

Uij, эВ	Q, безразмерный	Aij, эВ	Bij, 1/A	Cij, эВ*A6
Анион-Анион	0,68623	50211,7383	-5,5200148	74,7961121
Катион-Анион	0,68623	873,106567	-2,7838547	0
Катион-Катион	0,68623	0	0	0

Для сравнения методики и результатов настоящей работы мы выбрали эмпирические ПП для модели «точечных ионов» из работ других авторов (включая самые современные), которые условно можно отнести к разным категориям по методам их восстановления и использовавшимся опорным экспериментальным данным:

1. Аналитические расчеты в гармоническом приближении по упругим свойствам при близких к абсолютному нулю температурах (Catlow-81 [41], Busker-99 [42]).

2. Молекулярная динамика в ПГУ и низкотемпературные (до 1700К) данные по изотермической сжимаемости (Yamada-01 [50], Basak-03 [51]).

3. Расчеты методом статики решетки по энергиям образования и мигра­ции дефектов (Morelon-03 [43]).

На рис. 8.3 и рис. 8.4 приведено сравнение, температурных зависи­мос­тей периода решетки (полученных NPT МД-моделированием), а в табл. 8.4 – срав­нение упругих свойств (рассчитанных в гармоническом прибли­же­нии) для ПП из нашей работы, ПП из работ вышеуказанных авторов и экспери­ментальные данные [37][56]–[58].

Таблица 8.4.

Сравнение упругих свойств для всех ПП [41]-[43][50]-[51]

Источник \ Свойства	Eсвязи, эВ	K, ГПа	C11, ГПа	C12/C44, ГПа
Эксп-05, обзор [56]	-	166-205	194-228	76-77
Настоящая работа	49.5	123	217	76
Morelon-03 [43]	65.9	125	217	79
Basak-03 [51]	43.4	179	407	65
Yamada-01 [50]	45.7	181	418	62
Busker-99 [42]	104.5	259	532	122
Catlow-81 [41]	104.1	226	471	103
Эксп-76, обзор [57]	-	179-207	389-396	120/60
Оценка-63 [58]	106.7	-	-	-

Можно отметить, что потенциалы 1-ого типа [41] и [42], опирающиеся на старые экспериментальные данные [57][58] и наиболее грубое приближение «гармонических колебаний», демонстрируют наихудшее воспроизведение современных экспериментальных данных – как в области упругих свойств, так и по тепловому расширению.

Рис. 8.3. Сравнение температурных зависимостей периода решетки UO2, полученных NPT моделированием для всех ПП с экспериментом.

Рис. 8.4. Сравнение отклонений периода решетки для всех ПП от экспериментального.

Потенциалы 2-го типа, восстановлены в работах [50] и [51] само­согласованным методом молекулярно-динамического моделирования, что позволило напрямую учитывать тепловые (кинетические) свойства системы. Однако в качестве опорных экспериментальных данных в этих работах выбраны низкотемпературные и относительно плохо исследованные данные по изотермической сжимаемости (обе работы ссылаются на единственный экспериментальный источник [59]). Из сравнения видно, что эти ПП неплохо воспроизводят упругие свойства лишь из старых обзоров и тепловое расширение до 2000 К.

В работе Морелона [43] применен «противоположный» подход – ме­нее совершенный, но и менее ресурсоемкий метод статики решетки, не поз­во­ляющий учесть тепловые свойства системы, совмещен с опорными дан­ными по энергиям дефектов, не зависящими от температуры. Однако в оценках по энергиям дефектов имеется существенный разброс, например, для дефекта Френкеля рассматривается диапазон 3.5~5.4 эВ [37]. Соответствен­но, выбор этих данных в качестве опорных сопряжен с определенным риском. Тем не менее, можно отметить, что энергии дефектов, выбранные Морелоном, позволили получить потенциал, неплохо воспроизводящий современные экспериментальные данные по тепловому расширению и упругим свойствам, что можно считать косвенным подтверждением корректности выбранной Морелоном оценки для энергий дефектов, предложенной в работах Матцке [60]. К сожалению, в работе [43] параметры ПП приведены не полностью – потенциал анион-анион там кусочно-заданная функция, но в статье приведены точные параметры только для крайних отрезков, а детали их интерполяции отсутствуют. Поэтому корректно воспроизвести результаты работы [43] не представляется возможным, так что все данные для сравнения взяты из самой статьи.

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

Как видно из рис. 8.4, по отклонению модельного периода решетки от экспериментального удается воспроизвести тепловое расширение с точ­ностью 0.01 Å в диапазоне 500~2900К и с точностью 0.005 Å в диапазоне 900~2800 К, значительно сократив отклонение от экспериментальных данных на высоких температурах по сравнению с предыдущими результатами. Из табл. 8.3 видно, что упругие свойства также воспроизвелись с хорошей точностью и только значение объемного модуля упругости занижено по сравнению с современными экспериментальными данными.

По поводу большого разброса значений энергии связи (энергия превра­щения 1 моль кристаллического диоксида урана в ионизированный газ) Морелоном [43] уже отмечалась ее недоступность для прямого изме­ре­ния в эксперименте. В работе [58] оценка энергии связи получена на основе энер­гий ионизации кислорода и урана, которые также получены не экспери­мен­тально, а путем экстраполяции. Нужно отметить, что величина этой энергии, как показали расчеты, практически полностью определяется зарядами ионов, поэтому ПП 1-го типа с целочисленными 100%-зарядами имеют наибольшее значение энергии связи (близкое к [58]), ПП 2-го типа с 60%-зарядами имеют наименьшее значение, а ПП Морелона с ~81%-заря­дами и наши ПП с ~69%-зарядами дают значение между этими грани­цами.

Для проведения дальнейших более ресурсоемких исследований транс­портных свойств (массо - и теплопереноса) в диоксиде урана по результатам предыдущего сравнения были выбраны 3 лучших ПП: потенциалы из работ Морелона [43] и Базака [51] и полученный нами.

В табл. 8.6 приведены экспериментальные данные и результаты МД-моделирования с нашим ПП и ПП Базака для энергий образования и переноса точечных дефектов, а также данные из работы Морелона [43], рассчитанные методом статики решетки (СР). Детальный анализ механизмов массопереноса в диоксиде урана проведен в работах [54][61], поэтому здесь следует заметить следующее:

· результаты для ПП в пределах погрешности совпадают с оценками Матцке, полученными на основе большого набора экспериментальных данных;

· результаты Морелона занижены на величину порядка 0.3~0.5 эВ (что, впрочем, может быть обусловлено особенностями расчета методом статики решетки);

· результаты, полученные на основе ПП Базака существенно завышены (примерно в 1.5 раза) и попадают только на самые верхние грани­цы оценок [37][60].

Подробный анализ последних экспериментальных данных и теорети­ческое обсуждение теплопереноса в диоксиде урана можно найти в совре­менном обзоре [37].

На рис. 8.5 приведена экспериментальная температурная зависимость энтальпии из этой работы и расчетные кривые, полученные нами МД-моде­лированием на основе исследуемых ПП. Данные приведены в диапазоне 1500~3100 К, низкие температуры опущены (все кривые ложатся там на экспериментальные данные) чтобы различия на высоких были более явными. Можно отметить, что все кривые недостаточно хорошо воспроизводят экспе­риментальную зависимость на высоких температурах, тем не менее, резуль­таты с нашим ПП все же лежат заметно ближе в диапазоне начиная с 2200 К (максимальное отклонение наблюдается в области плавления и составляет 10% против 15% у предшественников).

Рис. 8.5. Сравнение температурных зависимостей энтальпии UO2, полученных NPT моделированием для ПП [9], [11] и нашего с экспериментом [12].

Рис. 8.6. Сравнение изобарной теплоемкости с учетом электронных и катионных дефектов для ПП [9], [11] и нашего с экспериментом [12]

Таблица 8.5.

Сравнение энергий дефектов для ПП [9], [11] и полученного по данным теплового расширения решетки

Дефекты анионной подрешетки диоксида урана	Эксп. обзоры [37][60]	МД наш ПП	МД ПП из [51]	СР [43]
E, эВ	±∆E, эВ	E, эВ	±∆E, эВ	E, эВ	E, эВ
Энергия Образования	Дефект Френкеля	3.5	±0.5	3.75	±0.25	5.4	3.17
Энергия Миграции	Вакансионная	0.6	±0.1	0.65	±0.06	0.69	0.33
Междоуз. (прямая)	-	-	1.33	±0.19	2.16	1.37
Междоуз. (непрямая)	0.9	±0.1	0.92	±0.02	-	0.65
Энергия Активации Диффузии	Вакансионная	2.4	±0.3	2.48	±0.13	3.39	1.92
Междоуз. (прямая)	-	-	3.21	±0.07	4.85	2.96
Междоуз. (непрямая)	2.7	±0.4	2.80	±0.19	-	2.24

На рис. 8.6 приведены графики изобарной теплоемкости, экспери­ментальные и модельные, полученные конечно-разностным дифференциро­ванием температурной зависимости энтальпии с учетом вклада электронных и катионных дефектов. Там также приведена дополнительная эксперимен­тальная зависимость из более старой работы Ральфа [62], на которой явно выражен пик теплоемкости при суперионном переходе (~2650 K), однако, погрешности измерений вблизи пика, указанные в работе, составляют 0.02 кДж/моль/К, т. е. перекрывают отличие от современных данных. Дальнейший выход теплоемкости на постоянное значение, т. е. отсутствие вклада дефектов после плавления кислородной подрешетки также не подтверждается более поздними экспериментами [37][63].