Размещено на http://www.allbest.ru/

25

Федеральное агентство по образованию

ГОУ ВПО «Алтайская государственная педагогическая академия»

Кафедра теоретических основ информатики

Курсовая работа

Создание модели возникновения Солнечной системы из межзвездного газа на базе численного моделирования с учетом гравитационного взаимодействия частиц

Выполнили студентки

363 группы

Черетун Инна Александровна

Мачалина Екатерина Викторовна

Научный руководитель

Алтухов Юрий Александрович

Барнаул 2009

Содержание

Введение

1. Общая постановка задачи

1.1 Модель образования Солнечной системы

2. Состав среды протопланетного диска Солнца

2.1 Уравнение состояния среды протопланетного диска

3. Численная двумерная модель протопланетного газопылевого диска

3.1 Уравнения газовой динамики в форме законов сохранения для криволинейной системы координат

4. Анализ результатов исследований

4.1 Модель образования планетной системы Солнца

4.2 Модель движения системы материальных точек

5. Потенциал межмолекулярного взаимодействия

5.1 Численный алгоритм

5.2 Краевые условия

5.3 Программа молекулярной динамики

5.4 Измерение макроскопических величин

Заключение

Литература

Введение

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

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

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

Строго говоря, эволюция протопланетного диска Солнца в ее полной и адекватной физическим процессам постановке должна решаться в рамках общей проблемы образования Солнечной системы. В этом случае эволюция протопланетного диска логически следует из определенной фазы эволюции Солнечной системы. В такой постановке задача является чрезвычайно сложной, и в настоящее время только намечаются пути ее решения. Решение этой задачи в полной постановке возможно методами численного моделирования на основе полных физически адекватных многомерных численных моделей образования Солнечной системы с использованием современных вычислительных систем, например, таких, как вычислительная система МВС-1000 Института прикладной математики им. М.В. Келдыша РАН. При отсутствии полных моделей актуальным становится построение упрощенных аналитических и численных моделей для отдельных этапов сценария образования Солнечной системы, способных правильно описать основные физические процессы на соответствующем этапе.

В настоящих исследованиях была использована приближенная аналитическая модель, впервые предложенная в работах [12,13] при исследованиях атмосферы вращающегося коллапсара.

Численные расчеты проводились на основе методик и программных средств двумерного программного комплекса, разработанного в Институте прикладной математики им. М.В. Келдыша РАН [14].

До сих пор мы изучали динамику систем, состоящих только из нескольких частиц. Однако на самом деле многие системы, такие как газ, жидкости и твердые тела, состоят из большого числа взаимодействующих друг с другом частиц. В качестве иллюстрации рассмотрим две чашки кофе, сваренного в одинаковых условиях. В каждой чашке содержится примерно 10²³- 10²⁵ молекул, движение которых с хорошей точностью подчиняется законам классической физики. Хотя межмолекулярные силы порождают сложные траектории каждой молекулы, наблюдаемые свойства кофе в каждом сосуде неразличимы и их сравнительно легко описать. Известно, например, что температура кофе, если его оставить в чашке, достигает комнатной и с течением времени больше не меняется. Как связана температура кофе с траекториями отдельных молекул? Почему она не зависит от времени, даже если траектории отдельных молекул непрерывно меняются?

Этот пример с чашкой кофе ставит нас перед проблемой: как можно, исходя из известных межмолекулярных взаимодействий, понять наблюдаемое поведение сложной многочастичной системы? Самый очевидный подход - решить эту задачу в лоб, моделируя на компьютере саму задачу многих частиц. Можно себе представить, что на каком-нибудь суперкомпьютере будущего будут решаться микроскопические уравнения движения для 10²⁵ взаимодействующих между собой частиц. На самом деле этот подход, называемый методом молекулярной динамики, применили к «небольшим» системам многих частиц, насчитывающим обычно от нескольких сотен до нескольких тысяч частиц, и он уже много дал для понимания наблюдаемых свойств газов, жидкостей и твердых тел. Однако детальное знание траекторий 10⁴ или даже 10²⁵ частиц ничего не даст, если мы не знаем, какие именно вопросы требуют выяснения. Какие основные свойства и закономерности проявляют системы многих частиц? Какие параметры нужно использовать для описания таких систем? К подобным вопросам обращается статистическая механика, и многие ее представления нашли отражение в этой главе. Тем не менее единственное, что требуется для работы над этой главой, это умение численно решать уравнения Ньютона, чем мы уже занимались, и некоторое знакомство с кинетической теорией.

1. Общая постановка задачи

Процессы образования протопланетных дисков и соответствующих им планетных систем существенным образом зависят от процессов эволюции космической системы, в которой рассматриваются эти явления. Это относится и к образованию планетных тел в Солнечной системе. Например, известно, что в межзвездных облаках не образуются изолированные планетные тела, более того, в них не наблюдается рост частиц пыли более 10^-5-10^-4 см [1]. Предполагается, что в облаках межзвездного пространства существуют процессы, препятствующие росту пылевых частиц. В одной из гипотез таким процессом, который «стабилизирует» размер частиц, является столкновение облаков в межзвездном пространстве [1].

Таким образом, образование протопланетного диска Солнца нельзя рассматривать вне зависимости от процессов образования Солнца как звезды, т.е. от модели образования Солнечной системы.

1.1 Модель образования Солнечной системы

В общих чертах модель образования Солнечной системы была принята в следующем виде.

1. Солнце и его протопланетный диск образовались путем единого процесса гравитационного сжатия вращающейся протосолнечной газопылевой туманности (аналогично, как это было предсказано Лапласом) [2], стр. 18; [3], стр. 99.

2. Формирование Солнца как звезды произошло за промежуток времени, равный примерно 0,1•10⁶ лет [3], стр. 101. Солнце за этот период аккумулировало около 90% своей массы. В это же время (одновременно с формированием Солнца) происходило образование протопланетного диска Солнца. На этой стадии Солнце окружено непрозрачной аккреционной оболочкой, которая поглощает интенсивное излучение молодого Солнца и переизлучает его в инфракрасном диапазоне.

3. Данные последних лет показывают, что коллапс межзвездной газопылевой туманности протекал таким образом, что, по крайней мере, часть этой туманности не была полностью испарена и гомогенизирована [4], стр. 26. На последующих этапах температура протопланетного диска Солнца падала и происходила конденсация первоначально высокотемпературного газа в той части, где ранее протекали процессы испарения.

4. Вторая стадия формирования Солнечной системы соответствует стадии Т Тельца до выхода Солнца на главную последовательность [3], стр. 100;

[5, 6], [2]. К началу второй стадии вокруг Солнца может сохраниться лишь незначительная по массе прозрачная часть аккреционной оболочки. Более значительная ее часть может находиться вдали от звезды в виде тора, окружающего звезду и входящего в состав протопланетного диска. На второй стадии идет более медленное формирование протопланетного диска Солнца, и эта стадия по ее продолжительности оценивается примерно в 10⁶ - 10⁷ лет [3], стр. 100; [5]; [2], стр. 207.

5. Солнечный ветер возникает на второй стадии. По разным источникам информации продолжительность солнечного ветра несколько различается [4, 3, 5], но, вероятно, ее можно оценить равной примерно 10⁶ лет.

6. Планетная система Земля-Луна образовалась из зоны протопланетного диска Солнца, находящейся на расстоянии около 1 а.е. от Солнца. Средние параметры среды этой зоны диска следующие: плотность ?10^-9 г/см³, температура ?400^оК [7], стр. 509.

2. Состав среды протопланетного диска Солнца

Для описания эволюции протопланетного диска Солнца весьма важен состав его среды.

По данным работ [2, 3] состав протопланетного диска Солнца на 98% состоит из газа, в котором обилия по массе молекулярного водорода, гелия и всех остальных веществ составляет соответственно 0,71; 0,28; 0,01. На пылевые частицы приходится по массе от 0,5 до 1,5%.

Одним из ключевых вопросов в эволюции протопланетного диска является поведение его пылевой компоненты, а именно: рост размеров частиц и возможность образования достаточно крупных тел, способных далее расти с помощью своего тяготения. Этот вопрос относится к числу наиболее сложных и не решенных до настоящего времени. По существу от решения этого вопроса зависит путь эволюции планетной системы Солнца. Если реализуется возможность независимого образования достаточно крупных твердых тел, дальнейший рост которых возможен за счет их тяготения, то это путь, который описывается моделью Шмидта-Сафронова [9], в противном случае может быть справедлива, например, капельная модель, предложенная Энеевым Т.М. и Козловым Н.Н. [8, 10, 11].

В межзвездных облаках размер частиц пыли не превышает 10^-4 - 10^-5 см [1], что обусловлено существованием процессов, которые ограничивают рост размеров частиц. Существуют ли такие процессы в протопланетном диске? Ответ на этот вопрос остается открытым. Ряд авторов утверждает, что в протопланетном диске Солнца может происходить рост размеров частиц при столкновении между собой за счет их слипания [2, 9]. Предлагаемые возможные механизмы слипания частиц пыли: ван-дер-ваальсовые силы; разные типы «радиационного» спекания [2], стр. 413; эффект холодной сварки [9], стр. 139 и другие. Произойдет ли слипание или дробление частиц при столкновении зависит от их относительной скорости. По данным работ [3, 7] в протопланетном диске Солнца частицы достигают распределения по размерам, в котором имеются и мелкие частицы размером около 1 мкм, поддерживающие высокую непрозрачность вещества диска, и крупные около 1 см. Средний размер пылевых частиц составляет несколько десятков микрометров.

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

Для описания поведения протопланетного диска необходимо знать уравнение состояния его среды.

2.1 Уравнение состояния среды протопланетного диска

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

Так, по данным работы [3], если на пылевые частицы приходится по массе около 1,5% вещества солнечного состава, то для такой среды молекулярный вес равен 2,53, а показатель адиабаты - 1,43. Описание протопланетного облака в приближении идеального газа дает достаточно точную картину поведения его некоторых средних характеристик, а именно тех, которые локально определяются газовой компонентой, даже в том случае, когда пылевая компонента конденсируется и превращается в твердое вещество.

3. Численная двумерная модель протопланетного газопылевого диска

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

3.1 Уравнения газовой динамики в форме законов сохранения для криволинейной системы координат

Используемая численная методика является вариантом методики для численного решения двумерных газодинамических течений в областях сложной формы с подвижными границами, разработанной Годуновым С.К. и Забродиным А.В. с соавторами [19]. Преимущества данного метода заключаются в том, что он позволяет выделить границы областей газа с разными физическими свойствами (например, контактные границы, ударные волны, ввести эйлерову границу и т.д.) и проследить за сложным движением этих границ. По положению границ, рассматриваемых как топологический «четырехугольник», строится два семейства координатных линий, разбивающие всю область течения на криволинейные четырехугольники - ячейки счетной сетки. Сетка изменяется во времени, подстраиваясь под характер движения границ. Поэтому ее конфигурация может быть достаточно сложной, а стороны ячеек существенно криволинейными.

Для построения конечно-разностного алгоритма запишем законы сохранения в интегральной форме в локальных криволинейных координатах.

4. Анализ результатов исследований

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

Как видно из результатов расчетов раздела 3, конфигурация протопланетного диска весьма чувствительна к зависимости угловой скорости вращения от цилиндрического радиуса (r). Следует отметить, что задавая массу диска (кольца) и функциональную зависимость Щ(r) могут быть получены как «плоские» протопланетные диски, так и тороидальные протопланетные кольца.

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

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

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

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

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

4.1 Модель образования планетной системы Солнца

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

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

2. В период окончания аккреционных процессов и перехода протопланетного диска в стадию формирования планетных тел в нем образуются крупномасштабные неоднородности в виде протопланетных колец, которые соответствуют зонам планет солнечной системы.

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

4. В процессе гравитационного взаимодействия эти сгущения сближаются, сталкиваются и образуют протопланетное облако или протопланетную систему планета и ее спутники. Эта стадия, по нашему мнению, хорошо описывается капельной моделью Энеева-Козлова.

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

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

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

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

Модель Энеева-Козлова образования планетной солнечной системы является альтернативной моделью по отношению к модели Шмидта-Сафронова. Отличие предлагаемой модели от модели Энеева-Козлова заключается в том, что протопланетные облака образуются не за счет слияния небольших газопылевых сгущений, образовавшихся в результате неустойчивости газопылевой среды протопланетного облака, а за счет фрагментации крупномасштабных неустойчивых образований в виде тороидальных протопланетных колец.

В рамках предлагаемой модели естественным образом может быть объяснено происхождение спутников планет. Основным положением при этом является то, что планеты и их спутники произошли из одного протопланетного кольца. Конкретный сценарий образования спутниковой системы может быть предложен в результате анализа экспериментальных фактов по данной системе - планета и ее спутники. Например, образование Луны в предлагаемой модели может произойти двумя основными путями. В одном случае на последнем этапе фрагментации протопланетного кольца образуются два протопланетных тела - протоземля и протолуна, находящиеся на близких орбитах. На заключительном этапе образования протосистемы Земля - Луна протоземля захватывает протолуну. В другом случае единое протопланетное облако, образовавшееся из протопланетного кольца, распадается на две части в результате его неустойчивости. Привлечение дополнительных экспериментальных фактов, например, таких как - Луна обеднена летучими элементами и другие, может дать возможность выбрать путь, по которому шло в действительности образование Земли и ее спутника Луны. Как видно, эти сценарии отличаются от столкновительной (giant-impact) модели, выдвинутой американскими учеными (Melosh and Sonet, 1986), по которой Луна образовалась в результате столкновения с Землей космического тела планетарных размеров.

4.2 Модель движения системы материальных точек

1. Задача. Имеется система N материальных точек с массами mi, i = 1, 2,..., N, взаимодействующих друг с другом с внутренними силами, на каждую из которых действует внешняя сила. Исходя из начальных координат xi, yi, и скоростей vxi, vyi, определите координаты и скорости материальных точек в последующие моменты времени.

2. Теория. Рассмотрим механическую систему из N материальных точек. Основной закон динамики:

где F'ij - внутренняя сила, действующая на i-ую материальную точку со стороны j -ой материальной точки, Fi - равнодействующая внешних сил, действующих на i - ую материальную точку со стороны тел не входящих в систему.

Дифференциальное уравнение второго порядка может быть представлено двумя дифференциальными уравнениями первого порядка. Имеем:

Зная внешние и внутренние силы, действующие на каждую материальную точку, можно определить их ускорения. Исходя из координат и скоростей точки в момент времени t, можно рассчитать координаты и скорости точки в следующий момент времени t + Д t.

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

3. Алгоритм.

1. Задают число материальных точек N, их массы mi, координаты xi, yi и проекции начальных скоростей vix, viy, силовое поле Fx = Fx (x,y), Fy = Fy (x,y), а также шаг по времени Д t.

2. Начало цикла по t. Дают приращение по времени: переменной t присваивают значение t + Д t.

3. Определяют проекции Fxi, Fyi равнодействующей всех внешних и внутренних сил, действующих на каждую i - ую материальную точку в момент t + Дt, и записывают их в массивы.

4. В цикле переобозначают координаты всех материальных точек, записывая их в массивы xx[i], yy[i].

5. В цикле перебирают все материальные точки и определяют проекции ускорения, скорости и координаты для каждой из них в момент t + Дt:

axi (t + Дt) = Fxi (t + Дt)/mi,

vix (t + Д t) = vix (t) + aix (t + Дt)Дt,

xi (t + Д t) = xi (t) + vix (t + Дt)Дt.

По аналогичным формулам вычисляют проекции на ось OY. Результаты записывают в массивы x[i], y[i], vx[i], vy[i].

6. Стирают изображения материальных точек в предыдущий момент времени t, координаты которых сохранены в массивах xx[i], yy[i].

7. На экране строят точки в следующий момент t + Дt, либо рисуют графики или выводят результат в числовом виде.

8. Возвращение к операции 2. Если цикл по t закончился, - выход из цикла.

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

program PROGRAMMA4;

uses dos, crt, graph;

const N=20; dt=0.01;

var m,Fx,Fy,x,y,vx,vy: array[1..N] of real;

Gd, Gm, i, j : integer;

ax, ay, F, l : real;

label Metka;

Procedure Sila; {Вычисление действующих сил}

label Metka;

begin

For i:=1 to N do

begin

Fx[i]:=0;

Fy[i]:=0;

end;

For i:=1 to N do

for j:=1 to N do

begin

if j=i

then

goto Metka;

l:=sqrt(sqr(x[i]-x[j])+sqr(y[i]-y[j]));

if l<2

then l:=2;

F:=-m[i]*m[j]/sqr(l);

Fx[i]:=Fx[i]+F*(x[i]-x[j])/l-(random(20)+10);

Fy[i]:=Fy[i]+F*(y[i]-y[j])/l+m[i]*10;

Metka:

end;

end;

Procedure Nach_uslov;

begin

Randomize; {Задание случайных координат и скоростей}

for i:=1 to N do

begin

m[i]:=2;

x[i]:=random(80)+60;

y[i]:=random(80)+60;

vy[i]:=random(30)+15;

vx[i]:=random(30)+15;

end;

end;

BEGIN

InitGraph(Gd, Gm, 'c:\bp\bgi'); {Инициализация графики}

Nach_uslov;

Repeat Sila;

for i:=1 to N do

begin

ax:=Fx[i]/m[i];

ay:=Fy[i]/m[i]; {Вычисление ускорений}

vx[i]:=vx[i]+ax*dt;

vy[i]:=vy[i]+ay*dt; {скоростей}

x[i]:=x[i]+vx[i]*dt;

if (x[i]<50)or(x[i]>350)

then vx[i]:=-vx[i];{отражение}

y[i]:=y[i]+vy[i]*dt; {координат}

if (y[i]<50)or(y[i]>350)

then vy[i]:=-vy[i];{от стенок}

for j:=1 to n do

begin

if (i<>j)and(abs(x[i]-x[j])<=2)and(abs(y[i]-y[j])<=2)

then

begin

x[i]:=x[i]+vx[i]*dt;

x[j]:=x[j]+vx[j]*dt;

y[i]:=y[i]+vy[i]*dt;

y[j]:=y[j]+vy[j]*dt;

vx[i]:=-vx[i];

vx[j]:=-vx[j];

vy[i]:=-vy[i];

vy[j]:=-vy[j];

end;

end;

end;

delay(5000);

setcolor(3);

Line(48, 46, 352, 46);

Line(48, 48, 48, 352);

Line(352, 46, 352, 352);

Line(352, 352, 48, 352);

setfillstyle(1,11);

bar(48,48,352,352);

setcolor(0);

for i:=1 to N do circle(round(x[i]),round(y[i]),2);

until KeyPressed;

CloseGraph;

END.

5. Потенциал межмолекулярного взаимодействия

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

U=V(r₁₂)+V(r₁₃)+…+V(r₂₃)+…== (5.1)

где V(r_ij) зависит только от абсолютной величины расстояния r_ij между частицами i и j. Парное взаимодействие вида (5.1) соответствует «простым» жидкостям, например жидкому аргону.

В принципе форму V(r) для электрически нейтральных атомов можно построить путем детального расчета, базирующегося на основных законах квантовой механики. Такой расчет очень труден и, кроме того, обычно бывает вполне достаточно в качестве V(r) выбрать простую феноменологическую формулу. Наиболее важными особенностями V(r) для простых жидкостей является сильное отталкивание для малых r и слабое притяжение на больших расстояниях. Отталкивание при малых r обусловлено правилом запрета. Иначе говоря, если электронные облака двух атомов перекрываются, некоторые электроны должны увеличить свою кинетическую энергию, чтобы находиться в различных квантовых состояниях. Суммарный эффект выражается в отталкивании между электронами, называемом отталкиванием кора. Слабое притяжение при больших r обусловлено главным образом взаимной поляризацией каждого атома; результирующая сила притяжения называется силой Ван-дер-Ваальса.

Одной из наиболее употребительных феноменологических формул для V(r) является потенциал Леннарда - Джонса

(5.2)

График потенциала Леннарда-Джонса показан на рис. 6.1. Хотя зависимость r^-6 в формуле (6.2) получена теоретически, зависимость r^-12

Рис. 5.1 График потенциала Леннарда-Джонса. Отметим, что r измеряется в единицах и V(r)-в единицах

Выбрана только из соображений удобства. Потенциал Леннарда-Джонса параметризуется «длиной» и «энергией» . Заметим, что при r= V(r) = 0. Параметр представляет собой глубину потенциальной ямы в точке минимума V(r); этот минимум расположен при расстоянии r=2^1/6. Заметим, что данный потенциал является «короткодействующим» и V(r) для r > 2.5 по существу равен нулю.

Удобно выражать длины, энергию и массу в единицах , и m, где m-масса частиц. Мы измеряем скорости в единицах (/m)^1/2, а время- в единицах . Для жидкого аргона параметры и потенциала Леннарда-Джонса составляют /k_B= 119.8 К и = 3.405 Е. Масса атома аргона равна 6.69*10^-23 г, и отсюда = 1.82*10^-12 с. Во избежание путаницы мы отмечаем все безразмерные или приведенные величины звездочкой. Например, приведенная двумерная плотность определяется соотношением р^*= р/².

5.1 Численный алгоритм

Теперь, когда у нас есть четко описанная модель системы многих частиц, необходимо познакомиться с каким-нибудь методом численного интегрирования для расчета траекторий каждой частицы. Уже в своих первых расчетах по моделированию уравнений движения Ньютона мы пришли к заключению, что устойчивость численного решения можно проконтролировать, следя за полной энергией и убеждаясь, что она не ушла от своего первоначального значения. Как можно было предполагать, алгоритмы Эйлера и Эйлера-Кромера не могут обеспечить сохранение энергии на временах, рассматриваемых при моделировании молекулярной динамики. К счастью, обычно нам не приходится привлекать сложный алгоритм. Интуитивно особенно привлекательным кажется алгоритм, определяемый формулами

x_n+1=x_n+V_nДt+a_n(Дt)², (5.3a)

V_n+1=V_n+(a_n+1+a_n) Дt. (5.3б)

Для упрощения обозначений мы записали этот алгоритм только для одной компоненты движения частицы. Алгоритм в виде (5.3) называется алгоритмом Верле в скоростной форме, и мы обсудим его в приложении 5А. В литературе по молекулярной динамике широко используется эквивалентная, но другая форма записи алгоритма Верле. Поскольку новая координата х_n+1 вычисляется с использованием не только скорости V_n, но и ускорения а_n, алгоритм Верле обладает «более высоким порядком» по Дt, чем алгоритмы Эйлера и Эйлера-Кромера. Новая координата используется для нахождения нового ускорения а_n+1, которое вместе с a_n используется для получения новой скорости V_n+1.

5.2 Краевые условия

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

Полезное моделирование должно включать в себя все характерные особенности рассматриваемой физической системы. Напомним, что конечной целью наших модельных расчетов является получение оценок поведения макроскопических систем, т.е. систем, содержащих порядка N=10²³-10²⁵ частиц. Рассмотрим сферический резервуар с водой. Доля молекул воды вблизи стенок пропорциональна отношению поверхности к объему (4рR²)/(4рR³/3). Поскольку N = с(4/3рR³), где с -плотность, доля частиц вблизи стенок пропорциональна N^2/3/N=N^-1/3, что при N?10²³ пренебрежимо мало. По сравнению с этим количество частиц которое можно изучать в моделях молекулярной динамики, составляет обычно 10² -10⁴, и доля частиц вблизи стенок не мала. В результате мы не можем провести моделирование макроскопической системы, помещая частицы в резервуар с жесткими стенками. Кроме того, если частица отражается от жесткой стенки, ее положение, а значит, и потенциальная энергия взаимодействия изменяются без какого-либо изменения в ее кинетической энергии. Отсюда присутствие жестких стенок означало бы, что полная энергия системы сохраняется.

Один из способов минимизировать поверхностные эффекты и более точно промоделировать свойства макроскопической системы заключается в использовании периодических краевых условий. Реализация периодических краевых условий для короткодействующих взаимодействий, таких как потенциал Леннарда-Джонса, хорошо знакома всем играющим в видеоигры. Рассмотрим сначала одномерный «ящик», содержащий N частиц, движение которых ограничено отрезком прямой длиной L. Концы отрезка прямой играют роль «стенок». Использование периодических краевых условий эквивалентно сворачиванию этого отрезка в кольцо (рис. 6.2). Заметим, что, поскольку расстоянию между частицами соответствует отрезок дуги, максимальное расстояние равно L/2.

В двумерном случае мы можем представить себе ящик, у которого противоположные ребра соединены так, что ящик превращается в поверхность тора (форма тороидальной камеры). Поэтому, если частица пересекает ребро ящика, она входит заново в противолежащее ребро. Заметим, что максимальное расстояние между частицами в х- и у-направлениях равно L/2, а не L.

Рис. 5.2

а - две частицы в точках х = О и х = З на отрезке длиной L = 4; расстояние между частицами равно З; б-отрезок преобразован в окружность; кратчайшее расстояние между обеими частицами на окружности равно 1.

То же самое можно представлять и по-другому, как проиллюстрировано на рис. 5.3. Предположим, что частицы 1 и 2 находятся в центральной клетке. Клетка окружена периодически повторяющимися собственными копиями.

Рис. 5.3 Пример периодических краевых условий в двумерном случае

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

Каждая копия клетки содержит обе частицы в тех же относительных положениях. Когда частица влетает в центральную клетку или вылетает из нее с одной стороны, это перемещение сопровождается одновременным вылетом или влетом копии этой частицы в соседнюю клетку с противоположной стороны. Вследствие использования периодических краевых условий частица 1 взаимодействует с частицей 2 в центральной клетке и со всеми периодическими копиями частицы 2. Однако для короткодействующих взаимодействий мы можем принять правило ближайшей частицы. Это правило означает, что частица 1 из центральной клетки взаимодействует только с ближайшим экземпляром частицы 2; взаимодействие полагается равным нулю, если расстояние до копии больше L/2 (см. рис. 5.3). Поскольку из данного правила следует, что мы можем визуально представлять себе центральную клетку в виде тора, это использование периодических краевых условий вместе с правилом ближайшей частицы было бы точнее назвать тороидальными краевыми условиями. Однако мы верны принятым традициям и называем эти условия периодическими краевыми условиями. Отметим, что использование периодических краевых условий означает, что все узлы ящика эквивалентны.

5.3 Программа молекулярной динамики

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

program md(input,output);

Begin

clrscr;

initial(x,y,vx,vy,N,nave,nset,Lx,Ly,dt,dt2);

pe:=0.0;

accel(x,y,ax,ay,N,Lx,Ly,pe);

time:=0.0;

pe:=0.0;

ke:=0.0;

xflux:=0.0;

yflux:=0.0;

virial:=0.0;

for iset:=1 to nset do

begin

for iave:=1 to nave do

Verlef(x,y,vx,vy,ax,ay,N,Lx,Ly,dt,dt2,virial,xflux,yflux,pe,ke);

results(N,nave,Lx,Ly,dt,virial,ke,pe,xflux,yflux,time);

end;

save_conf(x,y,vx,vy,N,Lx,Ly);

End.

Заметим, что х- и у-компоненты координат, скоростей и ускорений представлены массивами. Lх и Lу равны горизонтальной и вертикальной длинам прямоугольного резервуара. В переменных kе и ре накапливаются суммы для кинетической и потенциальной энергии; величины хflux, уflux и viriаl понадобятся для вычисления давления. Характер этих величин будет рассмотрен позднее.

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

Один из возможных вариантов начальных условий выглядит так: частицы помещают в узлах прямоугольной сетки и выбирают х- и у-компоненты скоростей случайным образом. В большинстве языков программирования имеется функция, которая генерирует последовательность «случайных чисел» на отрезке [0, 1]. Поскольку ЭВМ детерминирована, она не может вычислять последовательности действительно случайных чисел. Тем не менее компьютер может генерировать числа в совершенно неочевидном порядке, и, что касается нашей задачи, это различие не имеет значения. Разговор о последовательностях случайных чисел пойдет в гл. 11. В языке Тrue ВАSIС мы имеем возможность использовать функцию rnd, которая генерирует случайные числа в диапазоне 0? rnd<1. Случайные значения V_x на отрезке [-V_max,-V_min] формируются с помощью инструкции vx[i]:=random*vmax;

Количество частиц, линейные размеры системы и начальные координаты и скорости частиц задаются в подпрограмме initiаl. Заметим, что Lх и Lу измеряются в единицах у - параметра потенциала Леннарда - Джонса. Поскольку скорости выбираются случайно, их следует подправить, имея в виду, что начальный полный импульс в х- и у-направлении может просто не получиться равным нулю. Сетка в подпрограмме initiаl выбирается так, чтобы в начальный момент N = 12 частиц находились в левой половине ящика.

procedure initial(var x,y,vx,vy:component;

var N,nave,nset:integer;

var Lx,Ly,dt,dt2:real);

var

irow,icol,nrow,ncol,i:integer;

ax,ay,xscale,yscale,Mx,My,vscale,vmax,vxcum,vycum:real;

fname:string;

datain:text;

newconf:char;

begin

write('число частиц=');

readln(N);

write('размеры ящика Lx и Ly=');

readln(Lx,Ly);

write('шаг по времени=');

readln(dt);

dt2:=dt*dt;

write('число шагов по времени между усреднениями=');

readln(nave);

write('количество наборов усреднения=');

readln(nset);

write('новая конфигурация?(y/n)');

readln(newconf);

if(newconf='y') then

begin

write('число частиц в ряду=');

readln(nrow);

write('максимальная скорость=');

readln(vmax);

ncol:=N div nrow;

ay:=Ly/nrow;

ax:=Lx/ncol;

i:=0;

for icol:=1 to ncol do

for irow:=1 to nrow do

begin

i:=i+1;

y[i]:=ay*(irow-0.5);

if((irow mod 2)=0) then

x[i]:=ax*(icol-0.25)

else

x[i]:=ax*(icol-0.75);

vx[i]:=random*vmax;

vy[i]:=random*vmax;

vxcum:=vxcum+vx[i];

vycum:=vycum+vy[i];

end;

vxcum:=vxcum/N;

vycum:=vycum/N;

for i:=1 to N do

begin

vx[i]:=vx[i]-vxcum;

vy[i]:=vy[i]-vycum;

end;

end

else

begin

{write('имя файла с конфигурацией');

readln(fname);}

write('относительное изменение скорости=');

readln(vscale);

assign(datain,'novij.txt');

reset(datain);

readln(datain,N,Mx,My);

xscale:=Lx/Mx;

yscale:=Ly/My;

for i:=1 to N do

begin

readln(datain,x[i],y[i],vx[i],vy[i]);

x[i]:=x[i]*xscale;

y[i]:=y[i]*yscale;

vx[i]:=vx[i]*vscale;

vy[i]:=vy[i]*vscale;

end;

close(datain);

end;

end;

Рис. 5.4 Начальные условия, использованные в задаче 5.1 с параметрами, приведенными в подпрограмме initiаl

Затем мы реализуем алгоритм Верле для численного решения уравнений движения Ньютона. Обратите внимание на то, что скорость частично обновляется с использованием старого ускорения. Далее вызывается подпрограмма ассеl, которая вычисляет ускорение, используя новую координату, и скорость еще раз изменяется. Подпрограмма Vеrlеt обращается к подпрограмме реriоdiс, которая в свою очередь обеспечивает рассмотрение частиц только в центральной ячейке.

procedure Verlef(var x,y,vx,vy,ax,ay:component;

N:integer;

Lx,Ly,dt,dt2:real;

var virial,xflux,yflux,pe,ke:real);

var

i:integer;

xnew,ynew:real;

begin

for i:=1 to N do

begin

xnew:=x[i]+vx[i]*dt+0.5*ax[i]*dt2;

ynew:=y[i]+vy[i]*dt+0.5*ay[i]*dt2;

vx[i]:=vx[i]+0.5*ax[i]*dt;

vy[i]:=vy[i]+0.5*ay[i]*dt;

periodic(xnew,ynew,xflux,yflux,vx[i],vy[i],Lx,Ly);

x[i]:=xnew;

y[i]:=ynew;

end;

accel(x,y,ax,ay,N,Lx,Ly,pe);

for i:=1 to N do

begin

vx[i]:=vx[i]+0.5*ax[i]*dt;

vy[i]:=vy[i]+0.5*ay[i]*dt;

ke:=ke+0.5*(vx[i]*vx[i]+vy[i]*vy[i]);

virial:=virial+x[i]*ax[i]+y[i]*ay[i];

end;

end;

В подпрограмме ассеl для нахождения полной силы, действующей и каждую частицу, используется третий закон Ньютона. (Напомним, что нашей системе единиц масса частицы равна единице.) Обращение к подпрограмме sераrаtiоn обеспечивает, что расстояние между частицами никогда не будет больше Lх/2 в х-направлении и Lу/2 в у-направлении.

procedure accel(var x,y,ax,ay:component;

N:integer;

Lx,Ly:real;

var pe:real);

var

i,j:integer;

dx,dy,r,force,potential:real;

begin

for i:=1 to N do

begin

ax[i]:=0.0;

ay[i]:=0.0;

end;

for i:=1 to (N-1) do

for j:=(i+1) to N do

begin

dx:=x[i]-x[j];

dy:=y[i]-y[j];

separation(dx,dy,Lx,Ly);

r:=sqrt(dx*dx+dy*dy);

f(r,force,potential);

ax[i]:=ax[i]+force*dx;

ay[i]:=ay[i]+force*dy;

ax[j]:=ax[j]-force*dx;

ay[j]:=ay[j]-force*dy;

pe:=pe+potential;

end;

end;

Потенциал и сила вычисляются в подпрограмме f.

procedure f (var force,potential:real;

r:real);

var

ri,ri3,ri6,g:real;

begin

r:=1;

ri:=1.0/r;

ri3:=ri*ri*ri;

ri6:=ri3*ri3;

g:=24.0*ri*ri6*(2.0*ri6-1.0);

end;

Координаты частиц выдаются на экран в подпрограмме snapshot. Частота, с которой координаты частиц должны обновляться на экране, зависит от ?t.

procedure results (N,nave:integer;

Lx,Ly,dt:real;

var virial,ke,pe,xflux,yflux,time:real);

var

pflux,pvirial,E,T:real;

begin

if time=0.0 then

writeln('время,T,E,pflux,pviral');

time:=time+dt*nave;

ke:=ke/nave;

pe:=pe/nave;

E:=(pe+ke)/N;

t:=ke/N;

ke:=0.0;

pe:=0.0;

pflux:=((xflux/(2.0*lx))+(yflux/(2.0*Ly)))/(dt*nave);

xflux:=0.0;

yflux:=0.0;

pvirial:=(N*T)/(Lx*Ly)+0.5*virial/(nave*Lx*Ly);

virial:=0.0;

writeln(time:9:3,T:9:4,E:9:4,pflux:9:4,pvirial:9:4);

end;

5.4 Измерение макроскопических величин

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

Кинетическое определение температуры вытекает из теоремы о равнораспределении: каждый квадратичный член, входящий в выражение для энергии классической системы, находящейся в равновесии при температуре Т, имеет среднее значение , где k_B -постоянная Больцмана. Отсюда мы можем определить температуру Т системы в d-мерном пространстве соотношением

 (5.4)

где сумма берется по всем N частицам системы и d компонентам скорости. Скобки <...> обозначают усреднение по времени. Выражение (5.4) представляет собой пример связи макроскопической величины, в данном случае температуры, с временным средним по траекториям частиц. Заметим, что соотношение (6.4) справедливо только в том случае, если движение центра масс системы равно нулю-не хватает еще, чтобы движение центра масс резервуара изменяло температуру! В системе СИ температура Т измеряется в градусах Кельвина (К) и постоянная k_B= 1.38*10^-23 Дж/К. В дальнейшем мы будем измерять температуру в единицах е/k_B.

Еще одной тепловой характеристикой является теплоемкость при постоянном объеме

С_V = (дЕ/дТ)_V,

где Е - полная энергия. С_V есть мера количества тепла, необходимого, чтобы произвести изменение температуры. Поскольку теплоемкость зависит от размеров системы, удобно определить удельную теплоемкость на частицу, а именно с_V = С_V/N. Легче всего получить с_V путем нахождения средней потенциальной энергии при соседних температурах Т и Т + ?Т; с_V складывается из температурной зависимости потенциальной энергии и удельной теплоемкости, связанной с кинетической энергией, т.е. (d/2)k_B.

Чтобы определить среднее давление, предположим что частицы наxoдятся в резервуаре с жесткими стенками. Как известно, столкновения частиц со стенками резервуара приводят к тому, что на каждый элемент стенки действует средняя результирующая сила. Средняя сила, действующая на единичную площадку, и есть давление P газа. Давление можно найти, связывая силу в горизонтальном направлении со скоростью изменения соответствующей компоненты импульса частиц, ударяющихся о стенку.

Из тех же соображений можно получить среднее давление в отсутствие стенок. Поскольку давление в равновесном состоянии во всех направлениях одинаково, мы можем связать давление с общим переносом количества движения через элемент поверхности в любом месте системы. Рассмотрим элемент поверхности dА и предположим, что среднее количество движения, пересекающее в единицу времени нашу поверхность слева направо, равно S₊, а S_--среднее количество движения пересекающее поверхность в единицу времени справа налево. Тогда средняя сила F равна

F = S₊ - S_- (5.5)

и среднее давление определяется выражением

P = (5.6)

где F_n обозначает компоненту силы, нормальную к элементу поверхности. (В двумерном случае давление равно плотности потока импульса через единичный отрезок, а не через единичную площадку.) Что служит соответствующей мерой давления для потенциала Леннарда - Джонса?

Один из способов реализации этого метода определения среднего давления состоит в следующем. Представим себе, что на «ребрах» элементарной ячейки установлены четыре воображаемые поверхности. Тогда мы можем модифицировать подпрограмму реriоdiс и вычислять поток импульса за один шаг по времени.

procedure periodic(var xtemp,ytemp,xflux,yflux:real;

px,py,Lx,Ly:real);

begin

if xtemp < 0.0 then

begin

xtemp:=xtemp+Lx;

xflux:=xflux-px;

end;

if xtemp > Lx then

begin

xtemp:=xtemp+Lx;

xflux:=xflux-px;

end;

if ytemp < 0.0 then

begin

ytemp:=ytemp+Ly;

yflux:=yflux-py;

end;

if ytemp > Ly then

begin

ytemp:=ytemp-Ly;

yflux:=yflux+py;

end;

end;

Среднее давление можно найти, добавив в основную программу следующие инструкции:

xflux:=0.0;

yflux:=0.0;

Другой способ вычисления давлении вытекает из теоремы вириала:

, (5.7)

где r_i- координата i-й частицы, F_i-полная сила, действующая на частицу i со стороны всех остальных частиц, и сумма берется по всем N частицам. Вывод формулы (5.7) дан в приложении 6А. Чтобы с помощью вириала вычислить давление, добавим в подпрограмму Verlet следующую инструкцию:

virial:=virial+x[i]*ax[i]+y[i]*ay[i];

Для вычисления среднего давления добавим в подпрограмму оutput следующие инструкции:

Итак, мы встретили уже два примера-(6.4) и (6.7), содержащие связь макроскопических величин со средними по времени от функций координат и скоростей частиц системы. В равновесном состоянии эти средние не зависят от времени. В гл. 15 и 16 мы будем рассматривать второй вид усреднения - среднее по статистическим ансамблям. Связь этих двух методов усреднения кратко обсуждается в гл. 15.

При моделировании методом молекулярной динамики процесс установления равновесия зачастую может занимать существенную часть от общего времени счета. Как правило, в качестве начальных условий удобнее всего выбирать «равновесную» конфигурацию из какого-либо прежнего расчета, которая отвечает температуре и плотности, близким к требуемым. Следующая подпрограмма запоминает последнюю конфигурацию пропускаемого варианта и может быть использована в качестве последней подпрограммы, вызываемой программой md.

Следующие инструкции можно добавить в подпрограмму initiаl, с тем чтобы можно было воспользоваться какой-то прежней равновесной конфигурацией в качестве начальной конфигурации нового расчета.

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

ЗАДАЧА 1. Качественные свойства жидкости и газа

а. Во всех пунктах данной задачи для получения своей начальной конфигурации используйте следующие инструкции dаtа и rеаd. Инструкция RЕАD читает элементы данных из инструкций DАТА и присваивает их значения соответствующим переменным. Каждая инструкция READ читает очередной элемент из инструкции DАТА. Напишите короткую программу, в которой используются инструкции RЕАD и DАТА, и удостоверьтесь, что вы правильно понимаете работу этих инструкций.

Выберите N = 16, Lх = Lу = 6, rsса1е = 1 и ?t = 0.02. Чему равна приведенная плотность? Выберите nsnар = 20 и проведите моделирование на протяжении по крайней мере 200 шагов по времени. Чему равна начальная температура системы? Каков характер ваших снимков? Как долго система достигает равновесия? Каков ваш критерий равновесия? Вычислите полную энергию системы, температуру и давление. Не учитывайте при усреднении по времени неравновесные конфигурации.