ЗАДАЧА 16.17. Моделирование твердых дисков методом Монте-Карло

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

N = 16 и область размерами Lx = 4.41 и = 0.5*v 3*Lx. Вычислите приведенную плотность р* и сравните ее с максимально возможной плотностью. Чему равно среднее расстояние между частицами? Выберите начальные координаты частиц так, чтобы они образовывали треугольную решетку. Для начала имеет смысл величины пробных шагов взять равными dxmax = dymax = 0.1. Получается ли при таком выборе шагов коэффициент принятия равным примерно 50%? Сравните функции g(r) для р* = 0.95, 0.92, 0.88, 0.85, 0.80, 0.70, 0.60 и 0.30. Сохраняя неизменным отношение L^/L^ используйте последнюю конфигурацию предыдущего варианта в качестве начальной конфигурации нового варианта с более низкой р*. Подождите по крайней мере 20 шагов Монте-Карло на частицу, чтобы дать системе возможность прийти в равновесное состояние, и усредните функцию g(r) для nmcs - 100.

б.         Каково качественное поведение функции g(r) при высокой и низ-кой плотностях? Опишите, например, количество и высоту пиков,наблюдаемых у g(r).

в.         Используя соотношение (16.316), вычислите зависимость давле-ния от р*. Изобразите график отношения PV/Nk^r как функцию отр". «Объем» равен V = LL^. Является ли PV/Nk^ возрастающей илиубывающей функцией от р? Можно предположить, что при низкихплотностях система ведет себя как идеальный газ с точностью дозамены объема иа V - Mr. Сравните результаты, полученные для са-мой низкой плотности, с этим прогнозом.

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

"д. Вычислите «эффективный коэффициент диффузии» D путем опреде-леиия суммарного среднего квадрата смешения <R(t) > частиц после того, как равновесие достигнуто. «Время» t можно отождествить с числом шагов Монте-Карло на частицу. Находя D из соотношения D = = <R(t) >/4f, оцените D для значений плотности, рассмотренных в п. «а». Изобразите график зависимости p*D от р*. Как зависит D от р* для разреженного газа? Попробуйте определить область значений р*, где D резко падает. Наблюдаются ли какие-нибудь признаки фазового перехода?

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

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

Примем в качестве модели потенциал межчастичного взаимодействия Лениарда—Джонса

Свойства этого потенциала были описаны  в гл. 6.  Для атомов аргона параметры Леннарда—Джонса равны с = 3.405 Л и c/kB = 119.8 К. Удобно ввести безразмерные температуру и давление посредством соотношений Т = kBT/c и Р* = Ро^/е.

Применим теперь алгоритм Метрополиса к непрерывному потенциалу. Главное изменение, которое необходимо внести в программу hard_disk, состоит в замене подпрограммы overlap на подпрограмму test. В этой подпрограмме основной интерес представляет величина petest, которая представляет собой новую потенциальную энергию взаимодействия частицы с номером itrial, находящейся в пробной точке с координатами (xtrial, ytrial). Для простоты вычисление потенциальной энергии частицы itrial производится путем подсчета энергии взаимодействия со всеми остальными N -1 частицами, а не только с частицами в области действия потенциала (~2.3с). Величина petest сравнивается с peold — потенциальной энергией взаимодействия частицы itrial в точке (x(itrial), y(itrial)). Отметим, что величина ре должна вычисляться заранее и передаваться в подпрограмму test.

Подпись: