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

Молекулярно-динамические характеристики одномерной точно решаемой модели на различных масштабах

М.А. Гузев, Ю.Г. Израильский, М.А. Шепелов

Институт автоматики и процессов управления ДВО РАН, Владивосток, 690041, Россия

На основе точного решения для одномерной цепочки гармонических осцилляторов вычислены механические характеристики модели на различных масштабных уровнях. Получено представление для модуля Юнга как функции числа частиц на данном масштабе и стационарное распределение температуры.

Molecular dynamic characteristics of a 1D exactly solved model

at different scales

M.A. Guzev, Yu.G. Izrailskii, and M.A. Shepelov Institute for Automation and Control Processes FEB RAS, Vladivostok, 690041, Russia

Based on an exact solution for a 1D chain of harmonic oscillators, we have calculated mechanical characteristics of the model at different scale levels. Young’s modulus is represented as a function of the number of particles at a given scale and stationary temperature distribution is obtained.

1. Введение

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

Построение модели, отражающей выше указанную методологию, связано с преодолением противоречия

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

© Гузев М.А., Израильский Ю.Г, Шепелов М.А., 2006

процессов внутри нанообъектов. Это соответствует рассмотрению наноструктуры как гетероструктуры, то есть переходу к другому иерархическому масштабу.

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

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

2. Аналитическое решение

Уравнения движения для системы попарно взаимодействующих частиц имеют следующий вид [4]:

тх у = c(Xj+l - 2Xj + Xj-1), у = 1,..., п, (1)

где Ху — координата частицы j; т — масса частицы; с — постоянная взаимодействия. Если полагать х0 = иЬ, и — некоторая постоянная скорость, а координату хп+1 зафиксировать, то решение системы уравнений (1) соответствует простейшей дискретной модели одномерного кристалла, подвергаемого одноосному растяжению.

Пусть в начальном состоянии частицы находятся на расстоянии а друг от друга и их начальная скорость равна нулю. Выбираем начало отсчета совпадающим с начальным положением нулевой частицы. Тогда координата Ху представима в виде Ху = ]а + и у, где функция и у определяет смещение частицы j из положения рав-

новесия и является решением уравнения (1). При этом и у удовлетворяют условиям:

х0 иЬ, иу|ь=о 0, и у|ь =о 0 ип+1 °. (2)

Метод построения решения системы уравнений (1) представлен в [5]. Используя его, запишем решение для и у в виде:

иу =----- (п - } +1) -

-1 п +1

1 . nkj nk

— sinmkt sin-------------ctg---------

Юк n +1 2(n +1)

(3)

mk = 2X sin (nk/2(n +1)), X = y¡c/m.

Убедимся, что (3) удовлетворяет уравнению (1) при условии (2). Действительно, дифференцируя (3) по времени, имеем

n+

I ^

1 k=1

. nkj ^

юк sin mkt sin—ctg

пк

n +1 2(n +1)

(4)

С другой стороны, учтем, что

. nk (j +1) „ . nkj . nk (j -1)

sin———- - 2sin—— + sin-

n +1 n +1

. 2 nk . nkj

= -4sin ---------------sin——

n +1

n +1 n +1 Тогда отсюда и из (4) следует, что

Uj+1 - 2uj + Uj _х =

c n +1 k=1

. nkj nk

<x>k sm sin- ctg—--

n +1 2(n +1)

то есть справедливо уравнение (1). Очевидно, что

u.l = 0. с jlt =0 мени равно

uj | о = 0. Значение скорости в начальный момент вре-

jlt=0 n + 1

(n - j + 1) -

n +1 k=1

nkj sin-------ctg-

nk

n +1 2(n +1)

Воспользуемся формулой, представленной в [4], для нахождения суммы:

nkj nk

Е

k=1

sin-

-ctg

n +1 2(n +1)

=(n - j +1).

В результате получаем и1 о = 0.

Таким образом, решение уравнения (1) при условиях (2) представляет сумму монотонной и периодической функций, причем последняя является линейной комбинацией нормальных мод дискретной модели.

3. Пространственно-временное усреднение

Обозначим через L0 = а(п +1) длину кристалла в начальный момент времени. При изменении положения

нулевой частицы она равна Ь = Ь0 - иЬ. Относительное удлинение кристалла как целого определяется

є г = -

Ь — Ьо

—и

Ь0 а(п +1)

(5)

Локальное значение деформации є рассчитывается по формуле

є і =■

(6)

Величины є і являются дискретными характеристиками, а є Ь — характеристика континуальной механики, переход к которой определяется некоторой процедурой усреднения микрохарактеристик. Во-первых, это усреднение выполняется на некотором пространственном уровне. С этой целью выделим интервал А, содержащий N = — N1 +1 частиц, где частицы с номерами N1,

соответствуют началу и концу интервала А. Деформацию єА на масштабном уровне А определим как среднее значение локальной деформации на интервале А:

£є;

£д=^^ = —--------—. (7)

А N N

Во-вторых, следует учесть, что рассматриваемая задача имеет два временных масштаба. Первый гои( определяет время изменения расстояния между частицами при внешнем воздействии и он равен Ьои( = а/и. Второй масштаб ¿¡п( характеризует минимальное время изменения нормальной моды, которое по порядку величины равно Ь;п( ~ ^ш/е (молекулярная единица времени). В предельном переходе к континуальной модели величина а^е/ш совпадает со скоростью распространения о1 продольных возмущений (скорость звука) в материале, тогда Ь;п( ~ аи1 и параметр X = Ьои(/¿¡п( ~ иг/и >> 1 является большим. Если перейти к безразмерной переменной 5 = ь/ьоиА в аргументе нормальных мод, то имеем быструю переменную. По ней можно выполнить усреднение по времени, полагая

/ = Нш^]&/. (8)

тТ 0

Средние (/) зависят только от медленного времени, определяемого скоростью квазистатического нагружения. Из (3), (5), (7), (8), в частности, следует, что (еА ^ = е Ь. Предлагаемая процедура усреднения будет в дальнейшем использоваться также для вычисления значения напряжений. При этом получаемые характеристики являются естественными параметрами континуальной механики на интервалах времени, значительно превышающих молекулярную единицу времени.

4. Вычисление напряжения

Для расчета значений напряжения используется формула, приведенная авторами [3], полученная в рамках кинетической теории. В соответствии с предложенным

подходом сначала следует вычислить среднее значение напряжения стд, которое определяется частицами, попадающими в интервал А:

( 2 ^

°а=-У Ш + 2 X (руХу). (9)

2 i*]

А

Здесь Ху = х{ - Ху; Еу — сила, действующая между частицами с номерами i и j, вычисляемая через потенциал взаимодействия между ними и = и (Ху): ¥у = = -Эи/дХу, а суммирование выполняется по всем частицам, находящимся в интервале А.

Рассмотрим статический вклад в напряжение, определяемый вторым слагаемым, который для рассматриваемого потенциала равен

А і=N

X (хі+\ — хі)(хі+1 — хі —а)-

(10)

Обозначим

° і =—-77X

п + 1 к =1

1

------81П X

Юі-

X

. пк( і +1) . пкі

81П-----—------ — 81П-

п +1 п + 1

Подставляя (3) в (10), имеем

—С 2 =— Na X

А

пк

2(п +1)

єА+(єа) (2єА—(єа)) + Х(& і & і )

і = N1

(11)

Дальнейшее упрощение суммы в соотношении (11) основано на использовании усреднения по времени. Это приводит к тому, что ненулевой вклад в конечный результат дают гармоники с одинаковыми номерами, в результате получаем:

!,< ° і ° л=1 (п+т) )2 х

XX

к=1

1

Ю2

. пк(і +1) . пкі

2

31П

п +1

— 81П-

п +1

ctg

пк

2(п +1)

1

2А2

/ \ и

п +1

X

к=1

2 пк( і +1/2) 2 пк

^ ^——^2----

п +1

2(п +1)

Отсюда и из (11) сразу следует (^) = —£№2(( єа)+(єа)2) —

2А (п +1)2 к=1

2 пк( і +1/2) 2 пк

С082 ——-----------^2

(п +1)

2(п +1)

Формула (9) содержит кинетическое слагаемое (пропорциональное р 2/ш). При дальнейшем сопоставлении получаемых результатов с выводами континуальной

и і +1 — и і

— С

X

теории необходимо выделить конвективный вклад из кинетического члена. С этой целью находим (и^ = = и (п - у + 1)/(п +1) и

-Ь> - П+т І

nkj nk

cos rokt sin-------^ctg-

п +1 2(п +1)

Усредняя по времени, получаем с точностью до множителя безконвективный кинетический вклад в напря-

жение:

N.

m І (új-(új))2

Ni

m í w

N2 n

2| +1 І І І

2 I n + 11 j = N1 k=1

. 2 nkj 2 nk

sin —— ctg

п +1 2(п +1)

Объединяя все вычисления, имеем следующее выражение для напряжения:

N

{аь) = -^-г— ca({єд) + (єд) ) + N — 1

w2 1

+ -

2a(N — 1) (n +1)2

x

І

'=N x ctg

— cos2 %k(j + ^2) IX

j =N1 k =1 _

2 nk

П +1

П +1

2(n +1)

где учтено, что А = a(N -1). Чтобы выполнить суммирование по j, запишем

sin2 zkj - cos2 zk(j + V2) =

= -cos4p(co^-^pcos2zkj -sm^sin2zkj),

nk

Zk = —7 Й + 1

и воспользуемся справочными формулами

Ni ,, . cos(N + N2)zk .

X cos 2 zk j =-----:---------sin Nzk,

N

І

j=N1

sin zk

N

^ _ . sin(N1 + N2)zk .

І sin 2 zkj =-------:----------sin Nzk.

j=N1 sin zk

В результате получаем

(ал) = —TTT ca(( єл) + (єл) 2) +

N — 1

m 2N

+------------

2a( N — 1) S (N1, N 2, n) =

S (N1, N 2, n),

1

(12)

ХІ

k=1

(n +1)

sin Nzk

N sin( zk/ 2)

cos(N1 + N2 + V2)zkctg (zkl2)

5. Анализ результатов

Формула (12) определяет уравнение состояния одномерной системы, которая является континуальной моделью стержня длины Ь0. Уравнение состояния стержня, как известно [6], имеет следующий вид:

■[1 — a (T — То)] — Н,

2 = E

где 2 — напряжение в стержне; E — модуль Юнга; а — коэффициент линейного температурного расширения; T — температура. Пренебрегая эффектами второго порядка малости, перепишем формулу для 2 в виде:

L — La

+ а(Т — То)

(13)

Сравним формулы (12) и (13) в указанном приближении. С этой целью перейдем к статическому пределу в (12), полагая и ^ 0 и оставляя линейный вклад по деформации, тогда

, х Nеa , >

)-

С точки зрения физики статический предел соответствует Т = Т0, то есть (13) редуцируется к

2= E-

L—L

о

L0

Из последних двух соотношений сразу следует, что величина

E (N) =

Nca

N

N — 1 N — 1

имеет смысл модуля Юнга на масштабном уровне А. Дальнейшее сравнение (12), (13) дает

2ca2

■ S = а (Т — То).

(14)

Отсюда видно, что функция 51 характеризует распределение температуры вдоль стержня. Оно является неоднородным, поскольку зависит от положения точек N , N2, определяющих расположение и длину интервала усреднения А. Чтобы получить для 51 явное представление, выполним приближенное суммирование в (12). Примем следующие допущения

sin Nzk = Nzt.

sin — , cos-^ = 1,

2 2 2

'-к* 1у-к,

которые, без сомнения, справедливы для нижних частот (при малых К). Это приближение остается приемлемым и для больших частот, когда К близко к и, так как соответствующие слагаемые в 51 не могут существенно изменить величину суммы 51 из-за наличия множителя с^2( ЧІ2)> который мало отличается от нуля при к ~ п. Это позволяет записать 51 в виде

S = -

1

N(n +1) k=1

cos Zk(N1 + N 2 + V2)

8 n cos zk (N1 + N2 +1/2)

x

x

Полученный ряд сходится абсолютно. Распространяем суммирование до бесконечности (ошибка не превосходит величины O(1/и)), тогда

0 8 n cos kt ^ _ п(N1 + N2 +1/2)

S------------------------о > t _-2-•

п2 k_1 k n +1

Поскольку переменная t e (0, 2п), то, используя справочную формулу из [7] для суммы, получаем

S = -^-(3t2 -6m + 2п2).

3п

Формула (14) позволяет оценить неоднородность распределения температуры вдоль стержня. По порядку величины и/а-^с/т _ и_ в << 1, для металлов а ~ = 10-5 градуса, тогда (T-T0)~0.5-105SP, при этом |5| ~1.

6. Заключение

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

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

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

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

Работа выполнена при поддержке интеграционных проектов 06-П-СО-01-002, 06-П-УО-01-001 и частично при поддержке РФФИ, грант № 05-01-00618-а.

Литература

1. Панин В.Е. Методология физической мезомеханики как основа построения моделей в компьютерном конструировании материалов // Изв. вузов. Физика. - 1995. - Т. 38. - № 11. - C. 6-25.

2. Головнев И.Ф., Головнева Е.И., Конев А.А., Фомин В.М. Физическая

мезомеханика и молекулярно-динамическое моделирование // Физ. мезомех. - 1998. - Т. 1. - № 2. - С. 21-33.

3. Головнева Е.И., ГоловневИ.Ф., Фомин В.М. Особенности примене-

ния методов механики сплошных сред для описания наноструктур // Физ. мезомех. - 2005. - Т. 8. - № 5. - С. 47-54.

4. Кур-чанов П.Ф., Мышкис А.Д., Филимонов А.М. Колебания железнодорожного состава и теорема Кронекера // ПММ. - 1991. -Т. 55. - Вып. 6. - C. 989-995.

5. ЛурьеА.И. Операционное исчисление и его приложения к задачам

механики. - М.-Л.: ГИТТЛ, 1950. - 432 с.

6. Румер Ю.Б., Рыгвкин М.Ш. Термодинамика, статистическая физика

и кинетика. - Новосибирск: Сиб. унив. изд-во, 2001. - 608 с.

7. ПрудниковА.П., БрыгчковЮ.А., Маричев О.И. Интегралы и ряды. -

М.: Наука, 1981. - 800 с.

8. Заславский Г.М., Сагдеев Р.З. Введение в нелинейную физику: от маятника до турбулентности и хаоса. - М.: Наука, 1988. - 368 с.