Вычисление упругих характеристик кристаллических решеток графена и алмаза с применением многочастичных потенциалов

Введение

В настоящее время большое внимание уделяется исследованию упругих свойств кристаллической решетки графена (монослой графита).

В данной работе вычисление упругих модулей кристаллических решеток графена и алмаза ведется аналитически и с помощью компьютерного эксперимента. При вычислении используется метод молекулярной динамики (ММД) [1] с применением многочастичного потенциала Терсоффа [2], [3], Терсоффа-Бреннера [4], [5], а также потенциала Бреннера второго поколения [8].

• Аналитические расчеты: И. Е. Беринский
• Компьютерный эксперимент: В. А. Цаплин

Аналитическое вычисление упругих модулей

Будет скоро добавлено

Описание компьютерного эксперимента

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

[math]\underline{v} (t + \tau / 2) = \underline{v} (t - \tau / 2) + \underline{w} (t) \tau[/math]

[math]\underline{r} (t + \tau) = \underline{r} (t) + \underline{v} (t + \tau / 2) \tau,[/math]

где [math]\tau[/math] – шаг интегрирования. Ускорение [math]\underline{w}(t)[/math] вычисляется через приложенную к частице силу.

Для нахождения силы, приложенной к одной частице, рассматриваемая область разбивается на квадратные ячейки со стороной, не меньшей [math]R + D[/math]. Рассматриваются все частицы, находящиеся в одной ячейке, либо в соседних ячейках к ячейке, где находится частица, для которой вычисляется приложенная сила.

Колебания атомов гасятся введением коэффициента трения в уравнения движения. При этом достигается однородное деформированное состояние со значениями [math]\varepsilon_{11} = 0.001[/math] (при растяжении вдоль оси [math]OX_1[/math]), и [math]\varepsilon_{22} = 0.001[/math] (при растяжении вдоль оси [math]OX_2[/math]). При вычислении используется несколько ячеек периодичности и ставятся периодические граничные условия.

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

На втором этапе механические напряжения в решетке вычисляются по формулам [6]:

[math] {{\underline{\underline{\tau\hspace{-0.5mm}}}\hspace{0.5mm}}}_i = \frac{1}{2V} \sum_{\alpha} \underline{F}_{\alpha}^i \underline{A}_{\alpha}^i = \frac{1}{2V} \sum_{\alpha} \underline{F}_{\alpha}^i (\underline{r}_{\alpha}^i - \underline{r}_i), [/math]

иначе

[math] {{\underline{\underline{\tau\hspace{-0.5mm}}}\hspace{0.5mm}}}_i = \frac{1}{V} \sum_{\alpha} \frac{\partial \Pi^i}{\partial \underline{A}_{\alpha}^i} (\underline{r}_{\alpha}^i - \underline{r}_i), \quad \Pi^i = \frac{1}{2} \sum_{j (\neq i)} V_{ij} [/math]

Здесь [math]{{\underline{\underline{\tau\hspace{-0.5mm}}}\hspace{0.5mm}}}_i[/math] – тензор механических напряжений для частицы (атом углерода) с номером [math]i[/math]. При однородном поле деформации находится средний тензор напряжений [math]({{\underline{\underline{\tau\hspace{-0.5mm}}}\hspace{0.5mm}}}_i)[/math] по всем частицам. [math]V[/math] – объем ячейки периодичности (в двумерной постановке – площадь). [math]\underline{F}_{\alpha}^i[/math] – векторный коэффициент, равный [math]\displaystyle 2 \frac{\partial \Pi^i}{\partial \underline{A}_{\alpha}^i}[/math], где [math]\alpha[/math] – номер соседней частицы к частице с номером [math]i[/math]. [math]\underline{A}_{\alpha}^i[/math] – вектор относительного положения соседней частицы: [math]\underline{A}_{\alpha}^i = \underline{r}_{\alpha}^i - \underline{r}_i[/math], где [math]\underline{r}_i[/math] – радиус-вектор частицы с номером [math]i[/math], [math]\underline{r}_{\alpha}^i[/math] – радиус-вектор соседней частицы ([math]\alpha[/math]). [math]V_{ij}[/math] – энергия, приходящаяся на одну связь, см. потенциал Терсофа.

В двумерной постановке (графен) коэффициенты упругости вычисляются с помощью соотношений между компонентами напряжений и деформаций:

[math]\sigma_1 = C_{11} \varepsilon_{11} + C_{12} \varepsilon_{22},\quad[/math]

[math]\sigma_2 = C_{12} \varepsilon_{11} + C_{22} \varepsilon_{22},\quad[/math]

[math]\tau_{12} = 2 C_{44} \varepsilon_{12}[/math].

При этом упругие модули выражаются через эти коэффициенты по формулам:

[math] \nu = \frac{C_{12}}{C_{11}},\quad E = \frac{C_{11}^2 - C_{12}^2}{C_{11}},\quad K = \frac{1}{2}(C_{11} + C_{12}). [/math]

Здесь [math]\nu[/math] – коэффициент Пуассона, [math]E[/math] – модуль Юнга, [math]K[/math] – коэффициент объемного сжатия.

В результате компьютерного эксперимента для графена получены значения, приведенные в таблице 1. Количество цифр соответствует точности вычисления.

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

[math] \begin{array}{l} \sigma_1 = C_{11} \varepsilon_{11} + C_{12} \varepsilon_{22} + C_{12} \varepsilon_{33},\quad % \sigma_2 = C_{12} \varepsilon_{11} + C_{11} \varepsilon_{22} + C_{12} \varepsilon_{33},\\ % \sigma_3 = C_{12} \varepsilon_{11} + C_{12} \varepsilon_{22} + C_{11} \varepsilon_{33},\\ % \tau_{12} = 2 C_{44} \varepsilon_{12},\quad % \tau_{23} = 2 C_{44} \varepsilon_{23},\quad % \tau_{31} = 2 C_{44} \varepsilon_{31}. \end{array} [/math]

Упругие модули выражаются по формулам:

[math] \nu = \frac{C_{12}}{C_{11} + C_{12}},\quad % E = \frac{(C_{11} - C_{12}) (C_{11} + 2 C_{12})}{(C_{11} + C_{12})},\quad % K = \frac{1}{3}(C_{11} + 2 C_{12}). [/math]

Ячейка периодичности кристаллической решетки алмаза представляет собой куб с гранью [math]4 a_d / \sqrt{3}[/math], где [math]a_d[/math] – межатомное расстояние кристаллической решетки алмаза в равновесном состоянии.

На одну ячейку периодичности приходится 8 атомов углерода (рис. 2), координаты которых в осях, параллельных ребрам куба имеют вид:

[math]r_1 (0, 0, 0)[/math],	[math]r_2 (a_d / \sqrt{3}, a_d / \sqrt{3}, a_d / \sqrt{3})[/math],
[math]r_3 (0, 2 a_d / \sqrt{3}, 2 a_d / \sqrt{3})[/math],	[math]r_4 (a_d / \sqrt{3}, 3 a_d / \sqrt{3}, 3 a_d / \sqrt{3})[/math],
[math]r_5 (2 a_d / \sqrt{3}, 0, 2 a_d / \sqrt{3})[/math],	[math]r_6 (3 a_d / \sqrt{3}, a_d / \sqrt{3}, 3 a_d / \sqrt{3})[/math],
[math]r_7 (2 a_d / \sqrt{3}, 2 a_d / \sqrt{3}, 0)[/math],	[math]r_8 (3 a_d / \sqrt{3}, 3 a_d / \sqrt{3}, a_d / \sqrt{3})[/math].

Объем такой ячейки равен [math]64\,a_d^3\,/\,3\sqrt{3}[/math], т.е. на один атом приходится объем, равный [math]8\,a_d^3\,/\,3\sqrt{3}[/math]. Результаты вычислений представлены в таблице 2.

Квазистатическое вычисление упругих модулей

Правильная формулировка?

В результате компьютерного эксперимента для графена получены значения, приведенные в таблице 1. Количество цифр соответствует точности вычисления.

Таблица 1. Результаты вычислений для решетки графена
Потенциал	K, Н / м	Е, Н / м	nu	C_11 = C_22, Н / м	C_12, Н / м
Терсофф	175.4	406.2	-0.158	416.7	-66.0
Бреннер	200.3	235.8	0.411	283.8	116.7
Бреннер 2	201.0	242.7	0.396	287.9	114.1

Таблица 2: Результаты вычислений для решетки алмаза
Потенциал	K, ГН / м^2	Е, ГН / м^2	nu	G = C_44, ГН / м^2	C_11 = C_22, ГН / м^2	C_12, ГН / м^2
Терсофф	425	1056	0.087	642	1074	102
Бреннер	485	291	0.400	386	623	415
Бреннер 2	442	1049	0.104	721	1075	125

Динамическое вычисление упругих модулей

Для вычисления упругих модулей динамическим способом находится период гармонических колебаний кристаллической решетки c заданной длиной волны. Рассматриваются два вида колебаний: продольные и поперечные.

В результате упругие модули получили значения, приведенные в таблицах 3, 4.

Таблица 3: Результаты динамических вычислений для решетки графена
Потенциал	K, Н / м	Е, Н / м	nu	C_11 = C_22, Н / м	C_12, Н / м	G = C_44, Н / м
Терсофф	170.6	397.1	-0.164	408.1	-66.9	237.5
Бреннер	200.4	235.7	0.412	283.9	117.0	83.44
Бреннер 2	201.4	242.8	0.397	288.3	114.5	86.92

Таблица 4: Результаты динамических вычислений для решетки алмаза
Потенциал	C_11 = C_22, ГН / м^2	G = C_44, ГН / м^2
Терсофф	1018	610.7
Бреннер	590	361