Решение некоторых задач наномеханики в квазистатической постановке Текст научной статьи по специальности «Физика»

УДК: 539.1

А.С. Кравчук

Московский государственный университет приборостроения и информатики (Москва)

РЕШЕНИЕ НЕКОТОРЫХ ЗАДАЧ НАНОМЕХАНИКИ В КВАЗИСТАТИЧЕСКОЙ ПОСТАНОВКЕ

Abstract

A molecular quasi-static problems are considered. Solution of 1 D-problems are obtained. An effect of self-compression is detected. This effect consists of the existence of non zero internal forces and displacements in a chain without external loads. The reason of such effect is the long distance atom interaction. A non uniqueness of the solution is detected for a chain after destruction. The friction tractions are calculated for a long distance interaction of two different chains. Obtained solutions are compared with Tomlinson's solution and detected, that a long distance interaction leads to an asymmetries of the friction tractions curves.

В предыдущих работах [1-2] был развит итерационный метод решения задач наномеханики в квазистатической постановке, основанный на использовании метода Ньютона-Рафсона, и решен ряд одномерных задач. В настоящей работе решены некоторые новые задачи и проведен их анализ. Выявлен эффект самосжатия, состоящий в том, что при отсутствии внешних воздействий в конечной цепочке возникают внутренние перемещения и усилия; данный эффект обусловлен дальнодействием. Установлено, что решение задачи о равновесии цепочки после разрушения может быть неединственным. Вычислены силы трения при взаимодействии двух цепочек с учетом дальнодействия. Дано сопоставление данного решения с решением Томлинсона. Установлено, что учет дальнодействия приводит к асимметрии кривых трения.

1. Основные уравнения задач наностатики

Напомним постановку задачи о силовом взаимодействии двух типов атомов. Система уравнений зависит от вида кристаллической структуры контактирующих тел. Отметим, что используемые в настоящее время методы формирования геометрических моделей различных кристаллических структур конечных размеров описаны в [3]. Будем помечать атомы первого топа индексом «1», второго - «2». Для математической формулировки задачи об определении напряженно-деформированного состояния и его эволюции необходимо выбрать потенциалы взаимодействия П^, П2, П12: потенциал П дает силы взаимодействия между атомами тела «1», П2 - между атомами тела «2». Потенциал П12 определяет взаимодействие между атомами первого и второго тел. В молекулярной динамике движение каждого атома описывается уравнениями Ньютона [3], в молекулярной квазистатике - уравнениями Лагранжа. В неоднородных структурах необходимо вводить две группы сил. Силы Fap первой группы определяют воздействие

на данный атом других атомов этой же группы; силы Fap - воздействие со стороны атомов другой группы. Уравнения квазистатического равновесия имеют вид:

N N(3-<) ~

Z Fy+Z F„= 0, i = 1,2 (1.1)

y=1, y^i p=1

Пусть ri - радиус-вектор i -й частицы, ui - перемещения частиц относительно выбранного

I r r I

начального состояния. Введем текущее расстояние между частицами rij — | ri - Tj | и

r r r

заметим, что Fij — F (riJ ), где сила F может принадлежать и первой, и второй группе,

I

равно как и перемещения u . Зависимость Fij = F(гу ) является нелинейной, причем эта

нелинейность учитывает переход от монотонного роста модуля силы к монотонному убыванию при монотонном увеличении расстояния между атомами (в работе рассматриваются только парные взаимодействия). Конкретный вид зависимости

Г Г

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

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

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

2. Решения некоторых одномерных задач

Рассмотрим задачу, схематически представленную на рис. 1.

Пронумеруем атомы первой цепочки целыми числами «0, 1, 2,...,N1», второй -целыми числами « N +1, N + 2,...,N + N2». Обозначим через х (с соответствующим индексом) координаты начального положения атомов, через и - их перемещения вдоль общей оси Ох. Будем предполагать, что атом номер «0» первой цепочки закреплен, а атом номер «N + М2» второй цепочки перемещается из положения хм + м -итіп в

положение хх + х + итах с некоторым шагом h. Заметим, что в принципе можно задать силу, действующую на атом « N + N2» второй цепочки. В этом случае надо использовать

особые методы прохождения точки, в которой рост силы с ростом перемещений сменяется падением силы при дальнейшем росте перемещений. Аналогичная проблема возникает при задании перемещения атома номер N + N2 и применении метода Ньютона для

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

Задача ставится таким образом: зная текущее значение перемещения атома « N + ^2», найти перемещения атомов с номерами і = 1,2,...,Л^; j = N +1,...,N + Ы2 -1.

Разрешающая система уравнений получается приравниванием нулю правой части уравнений (2.1). Для рассматриваемой одномерной задачи эта система нелинейная; объединив все неизвестные перемещения в вектор и , запишем ее в виде:

П М (Гу ) = М2 еХР(-а (Гу - Ге )) - еХР(-2а (Гу - Ге )) ]

(1.2)

—о—о—о—о—о—о о—о—о—о-> 0 12 м] Vі м1 + м2

Рис.1. Расчетная схема одномерной задачи

У(М)=0. (2.1)

Здесь в уравнений, соответствующих атомам первой цепочки,

і

элемент Уі вектора У равен FiR - FiL , где

^ = 20^ £ |ехр[-а1(хр - х -г{е1} + ир -иі)] - ехр[-2а!(хр - хі - г(1)) + ир -М)]}+

Р=1+1 (2.2)

N1+М2 4 '

+ £ 1ЄХР[-а12 (х9 - хі - Г?4 + ид - М )] - ехР[-2а12 (хд - х, - ) + Ид - М )]>

^=N +1

^ = 2°>1а1 £ |еХР[-а1(хі - хр - ГР + М - ир )] - ЄХР[-2а1(хі - хр - Ге(1)) + иі - ир )]1 (23)

р=0

иі, ир - перемещения атомов с номерами «і, р » соответственно.

Аналогичным способом строятся уравнения равновесия атомов с номерами «і» во второй цепочке УІ = ^ - F1L, где

^ = 202а2 £ {еХр[-а2 (х? - хг - Г?2 + ид - и, )] - еХр[-2а2 (х? - хг - ГЄ2)) + ид - и, )]} (2.4)

д=і+1

^ = 21 |еХР[-а2 (х, - хд - + иі - и9 )] - еХР[-2а2 (х, - х9 - ГЄ2)) + М - ид )]}+

а=Ы1+\

N (2.5)

+ 2012а12 £ {еХР[-а12(х^ - хі - Г?2 + и9 - иг )] - еХР[-2а 12 (хд - хі - - М )]}

^ І е ^

д=0

(12) ,шп0, Л1) г (2) •

е ’ е

Величина г выражается через , ге1;; для используемого в вычислениях

12

потенциала Морса ге(12) = (ге(1) + ге(2))/2; так же подсчитываются постоянные 012а (см. [3]).

Будем предполагать, что атом с номером «0» первой цепочки закреплен, а атом с номером «N + Ы2» второй цепочки перемещается из положения хМі +^ -итіп в положение

хМі +N + итаХ с некоторым шагом Л. Тогда вектор неизвестных

И = (и1, и2,...,+ м2- і)Г ; индекс «Т » означает операцию транспонирования, т.е. и -

вектор-столбец. Во втором варианте задается сила, действующая на атом «N + Ы2»

второй цепочки, следовательно, и = (м1, И2,...,иМі+^ )Т .

Для решения системы (2.1) используем метод Ньютона-Рафсона:

к+1,і+1 к+1,і Гшт/ / ^ 1-1 л к+1,і к+1,і Г /(1-1 л к+1,і //-> /гд

u , = u , -[ат/ +и ^ , = u , -и]к+,г •Au , (2 б)

У к + 1,і+1 ґ- і 1ч г к +1

где И - (і + 1 )-я итерация искомого вектора перемещений и ;

Auк+1,г+1 = uk+1,г+1 - uк,г+1. Заметим, что нелинейность системы (2.1), обусловленная

свойствами потенциала Морса, такова, что матрица производных [ат / аu]-\+1і = [А]-+ і на

некотором шаге нагружения может оказаться вырожденной. Систематический способ решения данной проблемы указан в работе [4] (см. ниже).

Рассмотрим следующие задачи.

2.1. Деформация однородной цепочки с заданной на правом конце силой, равной нулю. Выберем в качестве нулевого приближения состояние, в котором расстояния между соседними атомами равны равновесным расстояниям для пары атомов. Пусть цепочка состоит из атомов алюминия, для которого D = 54,94, а = 1,3588, г = 2,866; система единиц измерений выбрана так, как это предложено в [3]. Вычисления были проведены для 20, 40 и 80 атомов в цепочке, количество итераций было фиксированым и равно 200, что приводило к погрешности приближенного решения, вычисляемого как евклидова норма разности двух последовательных приближений порядка 10-14 -10-15. Некоторые из полученных результатов представлены на рис. 2-3. На рис. 2 показано распределение внутренних усилий для 40 (сплошная кривая) и 80 атомов алюминия (пунктирная).

Ж Р К,1

1 . . .

16

Рис. 2. Распределение внутренних усилий при нулевой нагрузке

А икл

К,1 >

3 5 0 1(5 32 48 (54 80

Рис. 3. Распределение перемещений в той же цепочке.

На рис. 2-3 K, L - номера атомов в цепочках соответственно из 40 и 80 атомов.

Из представленных результатов видно, что, во-первых, усилия внутреннего сжатия в цепочке без нагрузки постоянны для атомов, находящихся во внутреннем массиве (вне «погранслоев»), и не зависят от количества атомов в цепочке. Во-вторых, перемещения

вне погранслоев распределены по линейному закону, причем угол наклона соответствующих прямых не зависит от количества атомов. Очевидно, что сжатие свободной цепочки увеличивает тангенс угла наклона касательной к кривой «усилия -перемещения» - жесткость цепочки («модуль Юнга»).

2.2. Разрушение однородной цепочки с заданным на правом конце монотонно растущим перемещением. В однородной цепочке из N +1 атомов возможны 3 вида разрушения: отрыв атома с номером «1» от атома номер «0», отрыв атома номер «N » от атома номер «N +1» и, наконец, одновременный отрыв атома с номером «1» от атома номер «0» и атома номер «N» от атома номер « N +1». Все три варианта разрушения получаются при надлежащем выборе начального состояния: отрыв атома «1» от атома «0» получается при задании некоторого дополнительного зазора («дефекта») между атомами « N » и « N +1», для получения отрыва атома « N » надо задать зазор между атомами «0» и «1». Соответствующие распределения перемещений для N +1 = 50 показаны на рис. 4; I - номер атома в цепочке. Видно, что и после разрушения имеет место краевой эффект: существуют погранслои и внутренняя линейно-упругая часть (ядро), в которой можно использовать решение задачи о растяжении линейно-упругого стержня.

/[V и,

1 1 1 1 1

1 дг^ 1

Рис. 4. Виды разрушения цепочки

Была сделана попытка объяснить наличие трех решений путем анализа характера стационарности потенциала на этих решениях. С этой целью был вычислен потенциал исследуемой системы как функция шага нагружения, полученные зависимости показаны на рис. 5 для следующих входных данных: a =1,3588, О =54,94, Ге =2,866, что соответствует атомам меди; число атомов равно 50; минимальное и максимальное значения перемещений последнего атома (справа) равны соответственно -1.5138316 (что соответствует самосжатию) и 25; количество итераций равно 300; количество шагов по внешним воздействиям равно 400 (сплошная кривая), 800 (пунктирная) и 1200 (штрих-пунктирная).

п

ы

-5560

-5800

і А / * / / /

1 11 / і / / / / 4

і / 1 / # і і/ / і / / / /

■ І і і і І І V У / / / /у /

2/--' /

0

номер шага

Рис.5. Зависимость потенциала от нагрузки

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

Решения для неоднородной цепочки опубликованы в [1-2].

3. Вычисление сил трения (задача Томлинсона - Френкеля - Конторовой = ТФК).

Работы [5-6] были пионерскими по проблеме вычисления сил трения, исходя из законов межатомного взаимодействия. Рассматривались две цепочки атомов, при этом воздействие одной из цепочек на другую заменялось гармонической силой. Рассматривались динамика и квазистатика такой системы. Модель Френкеля-Конторовой отличалась от модели Томлинсона учетом не только силового взаимодействия различных цепочек, но и силовым взаимодействием атомов одной цепочки. Впоследствии оказалось, что такие модели дают удовлетворительные результаты при расчете активного элемента атомно-силового микроскопа [7].

Применим алгоритм настоящей работы для решения задачи ТФК. Рассмотрим для простоты однородную цепочку, показанную на рис. 1. Будем предполагать, что атом номер «0» закреплен. Перемещения всех остальных атомов неизвестны и подлежат определению. Внешним воздействием теперь будет сила, распределенная по цепочке по закону (как и в [5]; индекс «1» означает, что в разложении Фурье силы берется только первое слагаемое):

2p (х + U)

a

]

(3.1)

где U - параметр нагружения; остальные параметры в формуле (3.1) постоянны. Таким образом, к силам, фигурирующим в уравнениях (1.1), действующим на атом номер « І », добавится сила

F

TO,i

f sin[-P(ire + u+ Uk)] , a

(3.-)

где Ге - равновесное расстояние для пары атомов, k - номер шага нагружения. Заметим, что в модели Томлинсона сила сопротивления атома его смещению от состояния

Оказалось, что даже такая простая модель позволяет получить важные качественные эффекты. В частности, можно доказать, что поведение силы трения (3.3) определяется отношением h = Ге / а периодов контактирующих решеток. Например, если h - рациональное число или же h - иррационально, но атомы распределены случайным образом с равной вероятностью относительно равновесного состояния (после сдвига положения к началу координат), то сила (3.3) будет равна нулю. Ненулевую силу трения можно получить, заменив выражение (3.1) отрезком ряда Фурье.

Более интересным является тот факт, что при I = 2р/г/ ka < 1 равновесное

состояние системы устойчиво и определяется только величиной и ; если же I > 1, то для каждого значения и существуют два состояния равновесия, одно из которых неустойчиво; для существования решения необходимо выполнение условия: 1 + I оов(2рХ; / а) > 0 (требование положительности второй производной потенциальной энергии). Эти утверждения вытекают из анализа уравнения равновесия /-го атома:

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

Результаты вычислений, представленные на рис. 6-7 получены для следующих входных данных: число атомов равно 20; ит)п = 0, итах = 5, а = 2ге; число шагов по нагрузке равно 100; амплитуда ./1=9. Ввиду сильной нелинейности задачи стабильные результаты получаются для числа итераций порядка нескольких сотен. Сплошные кривые соответствуют модели Морса, пунктирные - модели Томлинсона; k - номер шага нагружения.

равновесия линейно зависит от смещения ui, " i : Fr г- = -kut, k = const; = xt - x0, где

x0 = ire - начальное положение i -го атома, xi - текущее. Сила трения по определению определяется как

k(xi - x0 ) = -f; sin(2pxt / a) .

(3.4)

Назовем локальной текущей силой трения по Томлинсону FJo0 силу, действующую на одиночный атом с номером I на шаге нагружения k :

Flo = f sin[a(lr, + Uk)].

(3.5)

Рис. 6. Полные силы трения

Рис. 7. Локальные силы трения

Видно, что учет нелокального взаимодействия приводит к заметной асимметрии кривых, соответствующих модели Морса. Кроме того, максимумы сил трения в модели Морса заметно больше, нежели в модели Томлинсона, хотя характер поведения обеих кривых с определенной степенью точности один и тот же.

Работа выполнена при частичной финансовой поддержке грантов РФФИ №08-0100349 и №07-06-00269.

Библиографический список

1. Кравчук А.С. О моделях и решении задач механики наноконтакта/ А.С. Кравчук// Математическое моделирование систем и процессов: сб. науч. тр. - 2007. - №15. - С.123-141.

2. Goryacheva I. Numerical solution of nanomechanics problems. Development of the hybrid method// I. Goryacheva, A. Kravchuk, P. Neittaanmaki// Proceedings of VIII World Conference on Computational Mechanics (WCCM VIII) and V European Congress on Computational Methods in Applied Science and Engineering (ECCOMAS V). June 30 - July 5. 2008. Venice, Italy.

3. Рит М. Наноконструирование в науке и технике. Введение в мир нанорасчета /М.Рит. - Москва - Ижевск:RCD, 2005. - 159 с.

4. Коробейников С.Н. Применение метода конечных элементов к решению нелинейных задач по деформированию и потере устойчивости атомных решеток:

препринт /С.Н.Коробейников// Институт гидродинамики им. М.А. Лаврентьева СО РАН -Новосибирск, 1997. - 33 с.

5. Tomlinson G.A. A molecular theory of friction/ G.A. Tomlinson// Phil. Mag. Series. -1929. - Vol.7. - P. 935-939.

6. Френкель Ю.И. О теории пластических деформаций и двойникования/ Ю.И. Френкель, Т. Конторова //ЖЭТФ. - 1938. - Т. 8. - С.1340.

7. Bharat Bhushan. Nanotribology and nanomechanics/ Bharat Bhushan. - Berlin. -Springer. 2006. - 1148 p.

Получено 25.08.2008

6З