Высокомолекулярные соединения (серия А), 2021, T. 63, № 3, стр. 221-233

ВВЕДЕНИЕ

Атомы углерода способны создавать многочисленные структуры, среди которых в последнее время большое внимание исследователей привлекает моноатомный кристаллический слой – графен [1–7]. Данный наноматериал демонстрирует высокую электро- и теплопроводность, рекордные жесткость и прочность на разрыв, колоссальную удельную поверхность и другие свойства, обеспечивающие большое разнообразие его потенциальных применений. Замечательные свойства графена позволяют использовать его в качестве нанонаполнителя в полимерных композитах для улучшения их характеристик. Добавление графена в полимерные матрицы заметно модифицирует их механические, тепловые и электрические свойства [8–11].

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

ЦЕПНАЯ МОДЕЛЬ НАНОЛЕНТЫ В ПОЛИМЕРНОЙ МАТРИЦЕ

Графеновая нанолента представляет собой полосу постоянной ширины, вырезанную из плоского листа графена. Для упруго изотропного графена продольная и изгибная жесткости наноленты слабо зависят от ее ориентации. Для определенности рассмотрим наноленту, края которой имеют структуру зигзаг (рис. 1). Модель цепи, описывающая продольные и изгибные движения наноленты, подробно представлена в работах [16, 17].

Пусть нанолента лежит в плоскости xy вдоль оси x. В цепной модели наноленты предполагается синхронное движение всех атомов, расположенных на одной поперечной линии (линии, параллельной оси y). В этом случае линию атомов можно считать одним узлом молекулярной цепи, расположенной в плоскости xz (рис. 1).

Однослойная нанолента графена в 2D модели задается как молекулярная цепочка из N₁ звеньев, имеющая гамильтониан

где ${{{\mathbf{u}}}_{n}} = ({{x}_{n}},{{z}_{n}})$ – двухмерный вектор, задающий положение n-го узла цепи, ${{M}_{1}} = 12{{m}_{p}}$ – масса одного узла (масса атома углерода, m_p – масса протона).

Потенциал

описывает взаимодействие соседних звеньев, a = = 1.228 Å – продольный шаг, K = 405 Н/м – продольная жесткость цепи, r_n – расстояние между соседними звеньями ($a = {{r}_{c}}\cos (\pi {\text{/}}6)$, ${{r}_{c}} = 1.418$ Å – длина валентной связи C–C в графене).

Потенциал

описывает деформацию n-го угла цепи θ_n, энергия ${{\varepsilon }_{c}} = 3.5$ эВ задает изгибную жесткость цепи.

Потенциал ${{W}_{1}}({{{\mathbf{u}}}_{n}},{{{\mathbf{u}}}_{k}}) = {{W}_{1}}({{r}_{{n,k}}})$, ${{r}_{{n,k}}} = {\text{|}}{{{\mathbf{u}}}_{n}} - {{{\mathbf{u}}}_{k}}{\text{|}}$, отражает слабые невалентные взаимодействия узлов цепи n и k. Энергию взаимодействия узлов можно найти численно как сумму невалентных взаимодействий атома углерода из одной поперечной линии наноленты со всеми атомами другой поперечной линии, усредненную по поперечным смещениям первого атома.

Энергия невалентного взаимодействия атомов углерода нанолент выражается через потенциал Леннард-Джонса

с энергией взаимодействия ${{\varepsilon }_{0}} = 0.002757$ эВ, равновесным расстоянием ${{r}_{0}} = 3.8085$ Å [18]. При использовании этого потенциала зависимость энергии взаимодействия узлов цепи от расстояния между ними ${{W}_{1}}(r)$ с высокой точностью аппроксимируется потенциалом Леннард-Джонса $(m,n)$с энергией взаимодействия ${{\varepsilon }_{1}} = 0.00832$ эВ, равновесным расстоянием ${{r}_{1}} = 3.605$ Å, степенями n = 11, m = 5 (рис. 2).

Рассмотрим наноленту, находящуюся в полимерной матрице – в кристалле макромолекул полиэтилена (CH₂)_∞. Будем считать, что зигзагообразные цепочки ПЭ двигаются как твердое тело и всегда лежат поперек наноленты параллельно ее поверхности (рис. 1). Тогда каждую макромолекулу ПЭ можно заменить эффективным объединенным атомом, находящимся в плоскости xz.

Гамильтониан (1) описывает энергию наноленты, приходящуюся на ее продольную полосу ширины ${{h}_{0}} = 1.5{{r}_{c}} = 2.127$ Å. Макромолекула ПЭ имеет форму плоского зигзага с шагом ${{r}_{p}} = 1.53$ Å (длина валентной связи CH₂–CH₂) и углом α = 110°. Продольный шаг зигзага ${{a}_{y}} = {{r}_{p}}\sin (\alpha {\text{/}}2)$, ширина зигзага ${{a}_{x}} = {{r}_{p}}\cos (\alpha {\text{/}}2)$. Мы будем рассматривать только смещения зигзага как твердого тела. В этом случае на продольное расстояние h₀ приходится масса атомов зигзага ${{M}_{2}} = 14{{m}_{p}}{{h}_{0}}{\text{/}}{{a}_{y}}$ = = 23.76m_p.

В плоскости xz каждой макромолекуле ПЭ соответствует частица массы M₂. Чтобы определить взаимодействие между такими частицами, найдем зависимость энергии взаимодействия двух параллельных транс-зигзагов от расстояния между ними. При вычислении энергии взаимодействия воспользуемся приближением объединенных атомов, когда группы атомов CH₂ соответствуют частицам, взаимодействия которых задаются потенциалом Леннард-Джонса (4) с энергией взаимодействия ${{\varepsilon }_{0}} = 0.004052$ эВ и равновесным расстоянием ${{r}_{0}} = 4.5$ Å [19]. При вычислении энергия взаимодействия усреднялась по сдвигам зигзагов вдоль их продольных осей и по поворотам вокруг них.

Расчеты показали, что усредненная энергия взаимодействия транс-зигзагов ${{W}_{2}}(r)$ с хорошей точностью описывается потенциалом Леннард-Джонса (5) с показателями $n = 11$, $m = 5.5$. При нормировке энергии на длину цепи ${{h}_{0}}$ энергия взаимодействия ${{\varepsilon }_{1}} = 0.0325$ эВ, равновесное расстояние ${{r}_{1}} = 4.62$ Å (рис. 2).

Полимерная матрица, окружающая наноленту, моделируется двухмерной решеткой частиц массы ${{M}_{2}}$ с гамильтонианом

где $N = {{N}_{1}} + {{N}_{p}}$, N_p – число частиц (число молекул ПЭ). В 2D модели каждый узел решетки, имеющий координаты ${{{\mathbf{u}}}_{n}} = ({{x}_{n}},{{z}_{n}})$, соответствует положению проекции центра n-й зигзагообразной цепи ПЭ на плоскость xz (рис. 1).

Взаимодействие наноленты (узлов цепи) с полимерной матрицей описывается суммой

Вычисления показали, что потенциал взаимодействия узлов цепи с частицами полимерной матрицы ${{W}_{3}}({{{\mathbf{u}}}_{n}},{{{\mathbf{u}}}_{k}}) = {{W}_{3}}({{r}_{{nk}}})$ тоже с высокой точностью аппроксимируется потенциалом Леннард-Джонса (5) с показателями n = 11, m = 5.5, энергией ${{\varepsilon }_{1}} = 0.0162$ эВ и равновесным расстоянием ${{r}_{1}} = 4.125$ Å (рис. 2).

Схематически 2D модель системы нанолента + полимерная подложка представлена на рис. 3. Здесь цепочка маленьких серых частиц соответствует наноленте графена, а крупные черные частицы – макромолекулам ПЭ. Гамильтониан системы имеет вид

Данная модель легко может быть обобщена на случай многослойной наноленты. Гамильтониан системы, K-слойная нанолента + полимерная матрица, будет иметь вид

Здесь N_k – число узлов в k-й цепи (наноленте), номер звена ${{n}_{k}} = {{N}_{1}} + ... + {{N}_{{k - 1}}} + n$, число узлов в K нанолентах ${{N}_{c}} = {{N}_{1}} + ... + {{N}_{K}}$, N_p – число молекул полимерной матрицы, общее число узлов системы $N = {{N}_{c}} + {{N}_{p}}$.

Для нахождения стационарного состояния многослойной наноленты в полимерной матрице нужно численно решить задачу на минимум энергии взаимодействия

где $N = {{N}_{c}} + {{N}_{p}} = {{N}_{1}} + ... + {{N}_{K}} + {{N}_{p}}$ – общее число узлов системы, E – потенциальная энергия системы (если в формуле (9) положить все скорости ${{{\mathbf{\dot {u}}}}_{n}} \equiv 0$, то энергия E = H).

СТАЦИОНАРНЫЕ СОСТОЯНИЯ УГЛЕРОДНЫХ НАНОЛЕНТ В КРИСТАЛЛИЧЕСКОМ ПЭ

Найдем основные состояния углеродных нанолент в кристалле ПЭ. Вначале проанализируем зависимость энергии основного состояния наноленты от расстояния ее расположения до поверхности кристалла. Для этого рассмотрим однослойную наноленту длины $L = ({{N}_{1}} - 1)a = 12.65$ нм (цепочку из N₁ = 104 звеньев) и поместим ее на поверхность прямоугольника из ${{N}_{p}} = 82 \times 60 - 28$ = = 4892 узлов (размер 2D кристалла ПЭ 37.16 × × 23.27 нм²).

Из численного решения задачи на минимум энергии (10) следует, что наноленте энергетически более выгодно располагаться внутри кристалла ПЭ, чем на его поверхности. На рис. 3 представлены основные состояния наноленты при ее расположении на разных слоях молекул ПЭ вблизи поверхности кристалла. Пусть ${{E}_{i}}$ – энергия основного состояния системы нанолента $ \in $ кристалл ПЭ при расположении наноленты в $i$-м слое молекул от поверхности кристалла ($i = 1,2,...$). Определим нормированную на длину L разницу энергий $\Delta {{E}_{i}} = ({{E}_{i}} - {{E}_{1}}){\text{/}}{{N}_{1}}$. Моделирование показало, что $\Delta {{E}_{i}}$ монотонно уменьшается с ростом номера слоя i. Так, $\Delta {{E}_{i}} = - 0.0381, - {\kern 1pt} 0.0418$, –0.0425, –0.0427 эВ для i = 2, 3, 4, 6. При $i > 4$ энергетический выигрыш от погружения наноленты в кристалл практически не меняется (нанолента перестает чувствовать поверхность кристалла при ее помещении в четвертый слой молекул).

Чтобы оценить возможность расслоения многослойных нанолент в полимерной матрице, поместим две наноленты длины $L = 12.65$ нм (число звеньев ${{N}_{1}} = {{N}_{2}} = 104$) в центр кристалла из ${{N}_{p}} = 82 \times 60 - 56$ = 4864 молекул (размер 2D кристалла ПЭ $37.16 \times 23.27$ нм²) (рис. 4). Рассмотрим, как зависит энергия основного состояния системы ${{E}_{i}}$ от числа слоев молекул кристалла между нанолентами $i = 0,1,2,...$ (значение i = 0 соответствует двухслойной наноленте внутри кристалла). Определим удельную разницу энергий ΔE_i = = $({{E}_{i}} - {{E}_{0}}){\text{/}}{{N}_{1}}$. Решение задачи на минимум энергии (10) показало, что энергия системы две наноленты $ \in $ кристалл ПЭ монотонно уменьшается с увеличением числа слоев молекул ПЭ между нанолентами. Так, $\Delta {{E}_{i}} = - 0.0354, - {\kern 1pt} 0.0417, - {\kern 1pt} 0.0429$, – 0.0431 эВ для i = 1, 2, 4, 14. При $i > 4$ разница энергий практически не меняется при увеличении числа слоев (наноленты перестают взаимодействовать друг с другом, если между ними находятся более четырех слоев молекул полимерного кристалла).

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

ПРОДОЛЬНОЕ СЖАТИЕ НАНОЛЕНТ В ПОЛИМЕРНОЙ МАТРИЦЕ

Рассмотрим наноленту длины $L = ({{N}_{1}} - 1)a$, лежащую вдоль оси x, сверху и снизу которой находится прямоугольный фрагмент кристалла ПЭ такой же длины из 30 слоев молекул (ширина слоя H = 11.4 нм) (рис. 5). Найдем вначале основное состояние такой системы (прямоугольный фрагмент кристалла, вдоль центра которого расположена многослойная нанолента) $\{ (x_{n}^{0},z_{n}^{0})\} _{{n = 1}}^{N}$ как решение задачи (10). При $L = 24.4$ нм (N₁ = 200) покрывающая наноленту полимерная матрица будет состоять из ${{N}_{p}} = 55 \times 60$, при $L = 36.7$ нм (${{N}_{1}} = 300$) – из ${{N}_{p}} = 82 \times 60$, при $L = 49.0$ нм (${{N}_{1}} = 400$) – из ${{N}_{p}} = 109 \times 60$ молекул.

Затем продольно сожмем эту систему вдоль оси x, положив для всех атомов системы новые координаты ${{x}_{n}} = (1 - d)x_{n}^{0}$. Далее зафиксируем положение атомов из левого и правого края системы, а также вертикальные смещения атомов из нижнего и верхнего края системы. Снова решим задачу на минимум энергии системы (10) с этими граничными условиями. Тем самым мы получим стационарное состояние продольно сжатой наноленты в полимерной матрице при относительном сжатии $d \geqslant 0$.

Моделирование показало, что сжатие наноленты вместе с полимерной матрицей может происходить по двум сценариям: продольное однородное сжатие, при котором все атомы системы равномерно смещаются вдоль оси сжатия, или неоднородное сжатие, при котором в центре наноленты происходит ее зигзагообразный локализованный изгиб, а полимерная матрица подстраивается под этот изгиб (рис. 5).

Зависимости изменения энергии системы $\Delta E = E(d) - E(0)$ от относительного сжатия d для однослойных нанолент разной длины представлены на рис. 6. Как видно, при однородном продольном сжатии выполняется закон Гука (энергия системы растет пропорционально квадрату коэффициента сжатия ΔЕ ~ d²) и не выполняется при неоднородном сжатии. Моделирование показало, что однородное продольно сжатое состояние наноленты с полимерной матрицей устойчиво только при сжатии $d < {{d}_{0}} = 0.030$ (пороговое значение сжатия d₀ не зависит от длины наноленты). При более сильном сжатии устойчивым будет только неоднородное состояние с локализованным изгибом наноленты. Неоднородное сжатие энергетически более выгодно при $d > {{d}_{1}}$. Критическое значение сжатия d₁ зависит от длины наноленты, но всегда ${{d}_{1}} < {{d}_{0}}$. Для длины L = 24.4 нм (число звеньев цепи ${{N}_{1}} = 200$) ${{d}_{1}} = 0.0233$, для L = 36.7 нм (${{N}_{1}} = 300$) ${{d}_{1}} = 0.0188$, для L = 49.0 нм (${{N}_{1}} = 400$) ${{d}_{1}} = 0.0154$. Чем длиннее нанолента, тем при меньшем относительном сжатии неоднородный сценарий ее сжатия с локализованным изгибом становится энергетически более выгодным (d₁ → 0 при N₁ → ∞).

На рис. 7 представлены зависимости изменения энергии системы $\Delta E$ от относительного сжатия d для многослойных нанолент. Здесь тоже закон Гука выполняется только при однородном сжатии. Однородно сжатое состояние системы будет устойчивым при $d < {{d}_{0}}$. Пороговое значение сжатия d₀ зависит от числа слоев наноленты K: ${{d}_{0}} = 0.030$ при K = 1 и ${{d}_{0}} = 0.018$ при K = 2, 3. При длине $L = 36.7$ нм (${{N}_{1}} = 300$) состояние с неоднородным сжатием системы энергетически более выгодно при $d > {{d}_{1}}$. Критическое значение сжатия d₁ тоже зависит от числа слоев наноленты: ${{d}_{1}} = 0.0188$ при K = 1, ${{d}_{1}} = 0.0153$ при K = 2 и ${{d}_{1}} = 0.0124$ при K = 3. Таким образом, в полимерной матрице увеличение числа слоев наноленты приводит к понижению устойчивости однородного продольно сжатого состояния (критическое значение сжатия d₁ снижается при росте К).

Отметим, что если нанолента полностью находится в полимерной матрице (ее концы не выходят на поверхность полимерной матрицы), то сжатие матрицы практически не приводит к сжатию наноленты. Для иллюстрации поместим наноленту длины L = 24.4 нм (число звеньев цепи ${{N}_{1}} = 200$) в центр прямоугольной полимерной матрицы, состоящей из ${{N}_{p}} = 82 \times 61 - 54$ молекул (размер матрицы $36.7 \times 23.6$ нм²). Моделирование показывает, что при $d \leqslant 0.1$ продольное сжатие матрицы практически не приводит к сжатию наноленты. Так, при относительном сжатии матрицы d = 0.1 относительное сжатие наноленты ${{d}_{c}} = 0.003$, т.е. здесь всегда ${{d}_{c}} \ll d$.

ПРОДОЛЬНОЕ СЖАТИЕ НАНОЛЕНТ НА ПОВЕРХНОСТИ КРИСТАЛЛА ПЭ

Рассмотрим лежащий вдоль оси $x$ прямоугольный фрагмент кристалла ПЭ, состоящий из 61 слоя молекул (ширина слоя $H = 23.6$ нм). Поместим на нижнюю и верхнюю поверхности кристалла наноленты такой же длины, как сам фрагмент (рис. 8). Найдем вначале основное состояние такой системы (прямоугольный фрагмент кристалла, покрытый сверху и снизу нанолентами) $\{ (x_{n}^{0},z_{n}^{0})\} _{{n = 1}}^{N}$ как решение задачи (10). При длине покрывающих кристалл нанолент $L = 24.4$ нм (число узлов в цепи ${{N}_{1}} = 200$) кристалл будет состоять из ${{N}_{p}} = 55 \times 61$, при $L = 36.7$ нм (${{N}_{1}} = 300$) – из ${{N}_{p}} = 82 \times 61$, а при $L = 49.0$ нм (${{N}_{1}} = 400$) – из ${{N}_{p}} = 109 \times 61$ молекул.

Сожмем эту систему вдоль оси x, положив для всех атомов системы новые координаты x_n = $(1 - d)x_{n}^{0}$. Далее зафиксируем положение атомов из левого и правого края системы. Снова решим задачу на минимум энергии системы (10) с указанными граничными условиями. Мы получим стационарное состояние продольно сжатой системы с безразмерным относительным сжатием $d \geqslant 0$.

Моделирование показало, что сжатие нанолент вместе с полимерной подложкой тоже может происходить по двум сценариям: 1) продольное однородное сжатие нанолент и подложки, все атомы системы равномерно смещаются вдоль оси сжатия; 2) неоднородное сжатие нанолент, при котором в их центрах образуются локальные колоколообразные изгибы, направленные в противоположную сторону от подложки, а область под изгибом частично заполняется молекулами подложки (рис. 8). Первый сценарий реализуется при слабом сжатии, второй – при сильном.

Зависимости изменения энергии системы ΔE от ее относительного сжатия d для однослойных нанолент разной длины представлены на рис. 9. Как видно, закон Гука выполняется только при однородном сжатии. Моделирование показало, что однородное сжатие приводит к устойчивому состоянию наноленты на полимерной подложке только при $d < {{d}_{0}} = 0.017$ (величина порогового значения сжатия d₀ не зависит от длины наноленты). При более сильном сжатии (при $d > {{d}_{0}}$) устойчивыми будут только неоднородно сжатые состояния с локальными выпуклыми изгибами нанолент.

Неоднородное сжатие энергетически более выгодно при $d > {{d}_{1}}$. Критическое значение сжатия d₁ зависит от длины наноленты, но всегда ${{d}_{1}}$ < d₀. Для длины $L = 24.4$ нм (число звеньев цепи ${{N}_{1}} = 200$) ${{d}_{1}} = 0.0143$, для $L = 36.7$ нм (N₁ = 300) ${{d}_{1}}$ = 0.0136, для $L = 49.0$ нм (${{N}_{1}} = 400$) ${{d}_{1}} = 0.010$. Чем длиннее система, тем при меньшем ее относительном сжатии неоднородный сценарий сжатия с локальным изгибом нанолент становится более энергетически выгодным (d₁ → 0 при N₁ → ∞).

Зависимости изменения энергии системы $\Delta E$ от относительного сжатия $d$ для многослойных нанолент представлены на рис. 10. Здесь тоже закон Гука выполняется только при однородном сжатии. Однородное равномерно сжатое состояние системы будет устойчивым при $d < {{d}_{0}}$. Пороговое значение однородного сжатия ${{d}_{0}} = 0.017$ для однослойных нанолент, ${{d}_{0}} = 0.011$ для двухслойных и ${{d}_{0}} = 0.0085$ для трехслойных нанолент. Таким образом, при расположении нанолент на поверхностях кристалла ПЭ увеличение числа слоев нанолент приводит к понижению устойчивости продольного однородно сжатого состояния и способствует образованию локализованных выпуклых изгибов нанолент (d₀ уменьшается с ростом К).

Отметим, что ранее рассматривалось неоднородное сжатие однослойной наноленты, расположенной на плоской недеформированной (“жесткой”) подложке [20]. При сжатии нанолента на такой подложке может образовывать выпуклый изгиб (bubble) или складку – структуры, которые могут двигаться по наноленте, сохраняя свою форму. В рассмотренном нами случае поверхность кристалла ПЭ образует “мягкую” подложку, в области изгиба подложка меняет свою структуру, часть ее молекул частично заполняет пространство под нанолентой.

ЗАКЛЮЧЕНИЕ

Предложенная двухмерная модель нанолент графена в полимерной матрице позволяет промоделировать образцы полимера большого размера с включенными однослойными и многослойными нанолентами графена. С использованием этой модели показано, что нанолентам энергетически более выгодно находиться внутри полимерной матрицы (кристалла макромолекул полиэтилена), чем на его поверхности, а внутри матрицы многослойным нанолентам выгодно разделиться на отдельные слои. Результаты моделирования позволяют заключить, что вблизи температуры плавления наноленты с поверхности матрицы будут перемещаться внутрь, а внутри многослойные наноленты будут расслаиваться на однослойные.

Проведено моделирование сжатия нанолент графена, находящихся внутри полимерной матрицы. Независимо от числа слоев наноленты ее продольное сжатие вместе с матрицей происходит либо посредством однородного сжатия всей системы нанолента + полимерная матрица, либо путем неоднородного сжатия, при котором в центре наноленты возникает локальный зигзагообразный изгиб с подстройкой под него полимерной матрицы. Однородное сжатие ленты энергетически более выгодно при слабом сжатии матрицы, а неоднородное – при сильном. Сжатие наноленты, находящейся на поверхности матрицы, происходит по двум аналогичным сценариям независимо от числа слоев наноленты. В одном случае наблюдается однородное сжатие всей системы нанолента + подложка, а во втором – неоднородное сжатие, при котором у нанолент образуются локализованные колоколообразные изгибы, направленные в противоположную сторону от подложки, а области под изгибами частично заполняются молекулами подложки. Первый сценарий реализуется при слабом сжатии матрицы, второй – при сильном.

Вычислительные ресурсы предоставлены межведомственным суперкомпьютерным центром РАН.

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (код проекта 18-29-19135).