Журнал физической химии, 2022, T. 96, № 4, стр. 607-612

МОДЕЛИ И МЕТОДЫ ИССЛЕДОВАНИЯ

Для атомистического моделирования глинистых и цементных материалов методом классической молекулярной динамики (МД) в последние 15–20 лет широко используется силовое поле ClayFF [6, 8], которое обладает простотой и гибкостью в применении, но позволяет в то же время моделировать различные свойства материалов в хорошем согласии с имеющимися экспериментальными данными [5]. Новая модификация силового поля ClayFF (ClayFF-MOH) [5, 7] позволяет теперь более точно моделировать наночастицы конечного размера и связанные с ними поверхностные эффекты за счет параметризации в явном виде углов между атомом металла и гидроксилом (M–O–H) в структуре и на поверхности таких частиц. В структуре эттрингита присутствует два типа таких углов: Al–O–H и Ca–O–H.

Моделирование системы эттрингит – водный раствор солей NaCl/Na₂SO₄ проводилось в два этапа. На первом этапе на основе кристаллографических данных [1] была построена сверхъячейка эттрингита, состоящая из 4 × 4 × 2 элементарных ячеек в направлениях осей a, b и c кристалла соответственно. В итоге размер сверхъячейки составил 44.92 × 44.92 × 42.88 ų. Для этой модели рассчитывались структура и свойства собственно кристалла эттрингита: механические характеристики, количество и время жизни водородных связей, образованных всеми возможными парами донор-акцептор в кристалле. На втором этапе для моделирования систем эттрингит – водный раствор солей NaCl/Na₂SO₄ была построена аналогичная сверхъячейка размером 4 × 2 × 2, которую затем расщепили по кристаллографической плоскости (010) между столбиками эттрингита и привели в соприкосновение с отдельно приготовленным и уравновешенным слоем соответствующего водного раствора. Толщину этого слоя выбрали размером 30 Å, что позволило напрямую сравнивать полученные новые результаты с более ранними расчетами [8]. Таким образом, размеры симуляционных ячеек составляли 42.6 × 43.8 × × 55.6 ų для модели с водным раствором NaCl и 42.6 × 43.8 × 55.4 ų для модели с водным раствором Na₂SO₄. В первом случае в растворе содержалось 16 ионов Na⁺ и 16 ионов Cl^–, во втором – 16 ионов Na⁺ и восемь ионов ${\text{SO}}_{4}^{{2 - }}$. Количество молекул H₂O в заданном объеме раствора выбиралось исходя из плотности воды при нормальных условиях ~1 г/см³. В итоге моляльность водного растворов NaCl составила ~0.5 M, а водного раствора Na₂SO₄ – ~0.4 M. Такие концентрации примерно соответствуют экспериментально измеренным составам растворов в порах цементных бетонов [3, 8].

Все расчеты методом МД проводились с помощью пакета LAMMPS [9] с применением периодических граничных условий, а для описания всех межатомных взаимодействий в системах принимались параметры силового поля ClayFF согласно работам [5, 8]. Классические уравнения движения численно интегрировались по алгоритму Верле [10] с временным шагом 1 фс. Радиус обрезания короткодействующих взаимодействий составлял 12.5 Å, а для учета дальнодействующих электростатических взаимодействий применялся метод Эвальда [10]. Все построенные модели сначала приводились к термодинамическому равновесию при нормальных условиях (T = 298 K, P = = 1 бар) с применением термобаростата Нозе–Гувера для статистического NPT-ансамбля [10], затем равновесные расчеты проводились в NVT-ансамбле. Длина МД-траекторий составляла 500 пс для каждого типа ансамблей при моделировании свойств кристаллов эттрингита и 1 нс при моделировании поверхностей раздела с растворами солей. Выбор всех указанных параметров расчетов является вполне стандартным при моделировании такого рода систем [5, 10] и давно доказал свою надежность и эффективность.

ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

Анализ результатов вычисленных МД-параметров элементарной ячейки и плотности эттрингита с 26 молекулами H₂O показал хорошее согласие как с экспериментальными данными, так и с более ранними расчетами. Так, экспериментальное значение плотности эттрингита находится в диапазоне от 1.77 до 1.78 г/см³ [11], а расчет методом теории функционала электронной плотности (DFT) показал значение плотности 1.747 г/см³ [11]. Плотность эттрингита в настоящей работе получилась равной 1.811 г/см³ и 1.814 г/см³ для оригинальной и модифицированной версий силового поля ClayFF соответственно. Величины векторов a и c в экспериментах составили значения от 11.229 до 11.240 Å для вектора a и от 21.468 до 21.478 Å для вектора c [11]. Расчетные значения (методом DFT) векторов a и c составили 11.289 и 21.617 Å [11] соответственно. Вычисленные в настоящей работе значения вектора a составили 11.163 и 11.144 Å для оригинальной и модифицированной версий силового поля ClayFF соответственно. Значения вектора c составили 21.623 и 21.673 Å для оригинальной и модифицированной версий ClayFF соответственно. Следует отметить, что результаты расчетов с обеими версиями силового поля ClayFF практически схожи, что указывает на отсутствие больших изменений в равновесной структуре эттрингита при явном учете углов Al–O–H и Ca–O–H модифицированной версией ClayFF.

Результаты вычисленных механических характеристик эттрингита (объемный модуль, K, модуль сдвига, G, модуль Юнга, E, и коэффициент Пуассона, ν) и их сравнение с экспериментальными данными и результатами расчетов методом DFT [12] приведены в табл. 1. Воспроизведение механических характеристик модифицированной версией поля ClayFF улучшилось по сравнению с оригинальной версией, что обусловлено более корректным описанием колебания углов связи между катионами металлов и структурными гидроксилами.

Такое улучшение воспроизведения механических характеристик хорошо согласуется с более корректным описанием структуры и динамики водородных связей в новой модели. В кристалле эттрингита существует устойчивая сетка водородных связей, которая и отвечает за стабильность кристалла в направлении вектора c элементарной ячейки [11]. В эттрингите представлено довольно широкое разнообразие донорно-акцепторных пар, которые создают водородные связи. В качестве доноров могут выступать как атомы водорода структурных гидроксильных групп, так и атомы водорода молекул воды, а в качестве акцепторов – атомы кислорода как сульфат-ионов, так и молекул обоих структурных типов H₂O, а также атомы кислорода структурных гидроксильных групп. В итоге было проанализировано 11 возможных донорно-акцепторных пар, между которыми образуются H-связи. Для определения водородной связи между донором и акцептором использовался наиболее распространенный геометрический критерий [13, 14].

Время жизни водородных связей оценивалось как интеграл по времени автокорреляционных функций водородных связей [15]. В итоге по результатам вычислений наиболее прочные водородные связи для обеих версий ClayFF возникают между молекулами H₂O и сульфат-ионами, что отлично согласуется с расчетами методом DFT [16]. Однако наши расчеты показывают, что водородные связи, возникающие между гидроксильными группами и молекулами H₂O также играют существенную роль. Так, время жизни водородных связей между атомами водорода гидроксилов и кислорода молекул воды (пары Hh⋅⋅⋅Ow и Hh⋅⋅⋅Ow*, где * – молекула H₂O в канале эттрингита) значительно увеличивается для модифицированной версии ClayFF по сравнению с оригинальной версией ClayFF (~10 пс вместо ~0.3 пс). На рис. 1 показаны временные автокорреляционные функции водородных связей для пары Hh⋅⋅⋅Ow. Кроме того, для пары Hh⋅⋅⋅Ow на 20% увеличивается и среднее число водородных связей, приходящихся на один атом H гидроксила. Такое повышение количества и прочности водородных связей между гидроксильными группами и молекулами H₂O приводит к формированию более устойчивой сетки водородных связей в кристаллогидрате, что, в свою очередь, приводит к лучшему воспроизведению механических характеристик эттрингита модифицированной версией ClayFF (таблица 1).

Рис. 1.

Автокорреляционные функции водородных связей для пары Hh⋅⋅⋅Ow.

Для анализа взаимодействия водных растворов NaCl и Na₂SO₄ с поверхностью эттрингита были рассчитаны усредненные по времени профили плотности растворенных ионов и молекул H₂O как функции расстояния соответствующих атомов раствора до поверхности кристалла. На рис. 2 представлены профили плотности для ионов Na⁺ и Cl^– в направлении Z, перпендикулярном плоскости поверхности эттрингита. Для ионов Na⁺ наблюдается изменение положения пиков профиля плотности в первом гидратированном слое над поверхностью (не далее 3.5 Å от поверхности) для модифицированной версии ClayFF, что указывает на изменение внутрисферной адсорбции ионов. Такое изменение положения пиков ионов Na⁺ поверхностью кристалла можно объяснить более устойчивой структурой последнего по причине увеличения времени жизни водородных связей между гидроксильными группами и молекулами H₂O. Кроме того, более реалистичный учет ограничений на изменения углов связи М–О–Н сильно уменьшает вероятность прямого доступа ионов Na⁺ непосредственно к отрицательно заряженным атомам кислорода поверхностных гидроксильных групп. В то же время, поскольку структурные молекулы H₂O, находящиеся на поверхности кристалла и взаимодействующие с водным раствором, имеют меньшую подвижность в силу большей связанности с гидроксилами эттрингита за счет более прочных водородных связей, они также создают более плотную положительно заряженную поверхность эттрингита, в случае применения модифицированной версии ClayFF.

Рис. 2.

Профили плотности (ρ) ионов Na⁺ и Cl^– водного раствора NaCl.

Во втором приповерхностном слое (от 3.5 до 7.0 Å от поверхности) положение пиков профилей плотности ионов Na⁺ для обеих версий ClayFF практически совпадает, но высота пика для модифицированной версии ClayFF несколько меньше и он расположен все-таки несколько дальше от поверхности, что указывает на изменение относительной доли внешнесферной адсорбции. Такое изменение, очевидно, связано с обсуждавшимися выше структурными изменениями в первом приповерхностном слое раствора за счет большего отталкивания ионов Na⁺ поверхностью эттрингита для модифицированной версии ClayFF. Вследствие создания такой более плотной положительно заряженной поверхности с применением новой версии ClayFF, наблюдается и более сильная адсорбция ионов Cl^– на поверхности. Это заметно по увеличению высоты пика профиля плотности Cl^– на рис. 2 и его расположению немного ближе к поверхности. При этом для обеих версий ClayFF наблюдается только внутрисферная адсорбция ионов Cl^–.

Анализ профилей плотности для атомов водного раствора Na₂SO₄ показал, что поведение ионов Na⁺ ничем не отличается от поведения ионов Na⁺ в водном растворе NaCl, обсуждаемом выше. Для ионов (SO₄)^2–, как и для ионов Cl^–, сохраняется положение пиков при некотором увеличении их высоты при использовании модифицированной версии ClayFF, что в свою очередь указывает на более сильную внутрисферную адсорбцию в этом случае. Анализ профилей плотности молекул H₂O кристаллической фазы и водных растворов также показал, что для модифицированной версии ClayFF наблюдается более устойчивый и плотный молекулярный слой раствора на поверхности эттрингита, что в итоге приводит к большей гидрофильности поверхности и изменению внутрисферной адсорбции молекул H₂O на поверхности.

Кроме одномерных профилей плотности атомов в межслоевых пространствах эттрингита, нами были проанализированы усредненные по времени двумерные распределения плотности каждого типа атомов на кристаллической поверхности (контурные карты). На рис. 3 показаны такие контурные карты атомов кислорода сульфат-ионов на поверхности эттрингита из расчета с водным раствором Na₂SO₄. Вертикальными пунктирными линиями указаны границы каналов эттрингита, у которых расстояния между пунктирными линиями наименьшие. В каналах эттрингита указаны кислороды ионов (SO₄)^2–, принадлежащие кристаллической фазе. Между каналами располагаются атомы кислорода ионов (SO₄)^2–, адсорбированных из водного раствора Na₂SO₄.

Рис. 3.

Плотность распределения атомов кислорода сульфат-ионов (сплошные жирные контуры) на поверхности эттрингита из расчета с водным раствором Na₂SO₄. Слева – для оригинальной версии ClayFF. Справа – для модифицированной версии ClayFF.

Из рисунков видно, что для модифицированной версии ClayFF наблюдается более стабильная структура кристаллической фазы, практически все кислороды ионов (SO₄)^2– располагаются в каналах эттрингита. Это одновременно указывает на большую устойчивость структуры эттрингита за счет формирования более прочной сетки водородных связей в самом кристалле. Как упоминалось выше, анализ одномерных профилей плотности ионов (SO₄)^2– водного раствора Na₂SO₄ показал, что наблюдается бóльшая адсорбция ионов (SO₄)^2– для модифицированной версии ClayFF. Контурные карты это также подтверждают: заметна бóльшая плотность атомов кислорода ионов (SO₄)^2– (рис. 3 справа). Стабильные положение и форма усредненных по времени контуров в то же время указывают, что ионы (SO₄)^2– водного раствора не диффундируют заметно по поверхности, а находятся в довольно локализованном состоянии на поверхности эттрингита. Анализ контурных карт поверхностной плотности других растворенных ионов и молекул H₂O растворов NaCl и Na₂SO₄ тоже подтверждает выводы, сделанные выше на основе одномерных профилей плотности частиц: при применении модифицированной версии ClayFF наблюдается бóльшая локализация молекул H₂O кристаллической фазы и изменяются доли внутрисферной и внешнесферной адсорбции ионов и молекул воды на поверхности.

Для исследования динамических свойств растворов рассчитывались двумерные коэффициенты диффузии молекул H₂O на поверхности эттрингита. Коэффициенты диффузии вычислялись через среднеквадратическое смещение атомов по соотношению Эйнштейна−Смолуховского [10, 17]. При этом поверхностный слой раствора условно делился на три зоны в перпендикулярном поверхности кристалла направлении: ближайшая приповерхностная зона раствора толщиной 3.5 Å; зона, следующая за приповерхностной, также толщиной 3.5 Å; “объемная” зона – оставшаяся часть раствора в нанопоре дальше 7 Å от поверхности. В случае раствора NaCl во всех трех зонах наблюдается небольшое снижение диффузионной динамики молекул H₂O для модифицированной версии ClayFF. Например, в приповерхностной зоне коэффициент двумерной диффузии для оригинальной версии ClayFF составил ${{D}_{{{{{\text{H}}}_{{\text{2}}}}{\text{O}}}}}$ = 1.15 × 10^–5 см²/с, а для модифицированной версии ClayFF – ${{D}_{{{{{\text{H}}}_{{\text{2}}}}{\text{O}}}}}$ = 9.64 × 10^–6 см²/с. Эти результаты качественно подтверждают сделанные выше выводы о бóльшей структурной стабильности поверхности эттрингита и усилении адсорбции ионов Cl^– и молекул H₂O этой поверхностью. В случае раствора Na₂SO₄ для первых двух зон, наоборот, наблюдается небольшое увеличение коэффициента диффузии молекул H₂O для модифицированной версии ClayFF. В частности, в первой зоне коэффициент двумерной диффузии для оригинальной версии ClayFF составил ${{D}_{{{{{\text{H}}}_{{\text{2}}}}{\text{O}}}}}$ = 1.07 × 10^–5 см²/с, а для модифицированной версии ClayFF – ${{D}_{{{{{\text{H}}}_{{\text{2}}}}{\text{O}}}}}$ = 1.18 × 10^–5 см²/с. Это небольшое увеличение диффузионной подвижности можно объяснить бóльшим ионным радиусом (SO₄)^2– по сравнению с Cl^–. Как было сказано выше, ионы (SO₄)^2– сильнее адсорбируются поверхностью эттрингита в модифицированной версии ClayFF и при этом занимают несколько бóльший объем в приповерхностной зоне, что, в свою очередь, приводит к бóльшей подвижности молекул H₂O, неадсорбированных поверхностью эттрингита.

Таким образом, более точный учет углового взаимодействия между катионами металлов и гидроксилами в структуре эттрингита позволяет предсказать более высокую степень адсорбции ионов растворов и молекул H₂O на поверхности, а также влияет на количественное описание их подвижности. Увеличение адсорбции и изменение подвижности ионов и молекул растворов объясняется формированием более устойчивой и прочной сетки водородных связей как в кристаллической структуре эттрингита, так и на его поверхности. При этом явный учет углового взаимодействия M–O–H не влияет на качество описания термодинамически равновесной структуры кристалла, но позволяет лучше воспроизводить механические свойства последнего.

Работа выполнена в рамках Программы фундаментальных исследований НИУ ВШЭ в 2019–2021 годах. Исследование выполнено с использованием суперкомпьютерного комплекса НИУ ВШЭ [18].