Моделирование фазовых переходов первого рода методом молекулярной динамики тема автореферата и диссертации по физике, 01.04.14 ВАК РФ

РОССИЙСКАЯ АКАДЕМИЯ НАУК ОБЪЕДИНЕННЫЙ ИНСТИТУТ ВЫСОКИХ ТЕМПЕРАТУР

На правах рукописи

ЖАХОВСКИЙ Василий Викторович

МОДЕЛИРОВАНИЕ ФАЗОВЫХ ПЕРЕХОДОВ ПЕРВОГО РОДА МЕТОДОМ МОЛЕКУЛЯРНОЙ ДИНАМИКИ

01.04.14 - Теплофизика и молекулярная физика

Автореферат диссертации на соискание ученой степени кандидата физико-математических наук

Москва - 1996

Работа выполнена в Теоретическом отделе Объединенного Института высоких температур РАН.

Научный руководитель: чл.-корр. РАН, профессор,

доктор физико-математических наук Анисимов С,И.

Официальные оппоненты: доктор физико-математических наук

Багикиров А. Г. кандидат физико-математических наук Малышенко С. П.

Ведущая организация: Московский физико-технический институт

Защита состоится " " 1996 г. в часов на за-

седании Специализированного совета Д 002.53.03 при Объединенном Институте высоких температур РАН по адресу: 127412, Москва, Ижор-схая ул. 13/19

С диссертацией можно ознакомиться в библиотеке ОИВТ РАН. Автореферат разослан /^Л^АЛ 1995 г.

Ученый секретарь Специализированного совега, кандидат технических наук . А.Н.Давыдов

© Объединенный Институт высоких температур РАН, 1996

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность темы. Проблема испарения конденсированного вещества относится к числу классических. Она представляет принципиальный интерес и валена для многих приложений в научных исследованиях и технике. Интенсивному исследованию этой проблемы в последние годы способствовала разработка эффективных методов генерации высоких плотностей энергии, основанных на использовании лазерного излучения и мощных пучков частиц. Полученный таким образом обширный экспериментальный материал обычно интерпретируется на основе весьма упрощенных теоретических моделей (см., например [1, 2, 3]). При этом процессы, происходящие в газовой фазе, •где важную роль играют кинетические эффекты [4, 5, б], изучены гораздо более основательно, чем процессы в конденсированной фазе и межфазном переходном слое. Вместо анализа последних в большинстве работ используется простейшая одночастичная модель, рассматривающая испарени!' как выход наиболее быстрых атомов из потенциальной ямы, глубина которой равна средней энергии связи 110 • Чтобы вычислить поток испаренного вещества, нормальная к границе фаз компонента скорости у, интегрируется по "хвосту" максвелловского распределения с V, > \/2ио/т .

Такая модель является, очевидно, слишком упрощенной. В действительности, для атомов, находящихся в поверхностном слое и вносящих основной вклад в испаренный поток, энергия связи С/ не является фиксированной величиной, а зависит от ближайшего окружения данного атома. Она равна 1/0 только цо поряду величины. Существенно также следующее обстоятельство. Поверхностный атом с энергией связи II до перехода в газовую фазу совершает в среднем ехр II/кТ колебаний, т.е. остается в связанном состоянии в течение времени порядка г ехр и/кТ, где г по порядку величины равно обратной деба-евской частоте. При температурах кТ <<11 это время много больше характерного времени изменения энергия связи атома I/ (т.е. времени порядка т , за которое происходит заметная перестройка ближайше! о окружения поверхностного атома). Таким образом, процесс испарения должен рассматриваться как коллективный. Вследствие флуктуаций энергии связи в поверхностном слое образуются атомы клн групиы атомов, для отрыва которых необходима энергия на частицу, значительно меньшая, чем Ц0, к время жизни которых я связанной состоянии в поверхностном слое оказывается порядка т . Этя атомы с кинетической энергией порядка средней тепловой апосэт существенный

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

Целью работы является:

1. Создание высокоэффективного вычислительного алгоритма и программы моделирования двухфазных систем методом молекулярной динамики.

2. Численное моделирование термодинамически равновесной системы жидкость-пар я процесса испарения простой жидкости в вакуум.

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

4. Выяснение роли флуктуаций энергии связи атомов в поверхностном слое в процессе испарения и определение механизма испарения простых жидкостей. .

5. Определение вида кривых фазового равновесия (в частности кривой плавления) и зависимости коэффициентов переноса от температуры фазового перехода в системах с однородной потенциальной функцией.

Научная новизна работы состоит в следующем:

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

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

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

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

•

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

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

Апробация работы. Результат исследований были представлены на International Colloquium on Gas Kinetics in the Presence of Phase Transitions (Minsk, 1991), 3-rd. International Symposium on Subsecomi Thermophysics, (Yraz, Austria, 1902), нескольких семинарах Теоретического отдела ИВТ РАН.

Публикации. По материалам диссертации опубликовано о Печатных работ.

Структура и объем диссертации. Диссертация состоит нч внедения, четырех глав и приложения, содержит 121 страниц машинописного текста, 57 рисунков на 57 страницах, 9 таблиц, 52 наименований ис-

. пользованной литературы.

•» • -

СОДЕРЖ АНИЕ ДИССЕРТАЦИИ

Во введении обоснована актуальность исследуемой проблемы, сформулированы основные цепи работы, указаны оснонные моложе- , ния, ныносимые на защиту, а г<икже описана структур* днссерпицш.

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

Вторая глава посвящена разработке высокоэффективной методики моделирования методом МД двухфазных систем.

Рассматривается система из N частиц, помещенных в параллелепипед с размерами £х £ х2£,, называемый далее МД-ячейкой. Взаимодействие между г -ой и з -ой частицами описывается потенциалом Лениард-Джонса: <

и(гу)=4е[(а/г,7)12-(ег/гу)в]) (1)

где гу - расстояние между частицами, а е и о ~ параметры потенциала. Чтобы разделить систему на жидкую и газовую фазы по оси г , вводится дополнительный потенциал, действующий на систему со стороны одной из граней параллелепипеда - дна МД-ячейки г = — Ь:

щЫ=А£ь{(сг/(1 + 2^~(а/(Ь + г1))% (2)

Для моделирования испарения в вакуум я в замкнутый объем на противоположной грани ячейки г = Ь вводится отталкивательный потенциал, играющий роль крышки МД-ячейки

«с(г,-)= 46^/(1-^)]12 (3)

При закрытой крышке ес = е, при открытой - ес = 0. В плоскости ху на систему налагаются условия периодичности с периодом Ь.

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

шг( =^(гьг2,...,г/у) » = 1,2,..,^ (4)

где силы £ определяются дифференцированием потенциалов (1)-(3).

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

I «

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

хп+1=2 + . (5)

Здесь п- номер временного шага, к - шаг по времени, а величины а] приведены в таблице.

] 0 1 2 3 4 5 6 7

аг 88324 -121797 245598 —300227 236568 -117051 33190 -4)25

Чтобы упростить написание, в формуле (5) выписаны только ,г компоненты вектора г.

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

2л(*»+| + £(*■•+» +2*.)

Основное отличие (6) от обычных конечно-разностных оценок первой производной в использовании дня повышения точносчи вторых про изводных. Этот подход представляется естественным, иоскч-'ч.ку используются значения уже вычисленных для (5) ускорений. С.зедуе! отметить, что значения скоростей вычисленных по (б) на н-ом шше времени не используются при интегрировании уравнений движения на (п + 1) -ом шаги. Полому ошибка оценки (6) влияет только ни значения зависящих от скоростей физических величии (тина кинетической энергии) вычисленных ни п -ом шаге. Иными слонами, ошибка меюла (6)'.нигде не накапливается.

При выполнении расчетов удобно использовать безразмерные переменные. Общепринятым является выбор а в качестве едшишы ;ии ны, ей качество единицы энергии и величины"га 0\/т) 48г в качестве единицы времени. В дальнейшем, если не указаны явно размерности, будут нсчользопагься чти единицы. Например, температура и без размерных единицах будет выражаться форму лой: Т » 10

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

параметрами потенцила Л е! ш ард-Джон с а конкретного вещества. Например, для аргона можно йсцользоввать а = 3.405 • Ю-8 см, е = 1.653-Ю-14 эрг =119.8 °К; то*!да единица времени г = 3.114• Ю-13с.

Чтобы начать интегрирование уравнений (4), необходимо задать начальные координаты и скорости всех частиц. Для ¿того частицы рас- . полагались в узлах простой кубической решетки внутри;куба ЬхЬхЬ у дна МД-ячсЙжи. Скорости частиц выбирались равными по, модулю V — уГ/16 и случайно распределенными по направления)«, Заметим, что результаты расчета не зависят от способа задания начальных условий. Как показали проведенные МД-экспериментй, максвеллов-, ское распределение в двухфазной системе устанавливается гораздо быстрее, чем достигается равновесие между фазами, Чтобы привести систему в равновесное состояние с заданной температурой, на начальном этапе моделирования приходится применять масштабирование скоростей и демпфирование колебаний центра масс вдоль оси г . После этого начального этапа следует более короткий промежуточный этап, на котором система не подвергается внешним воздействиям и переходит в стационарное состояние. При этом происходит небольшое измененние температуры. Критерием завершения обоих этапов служит постоянство температуры вдоль оси г и малость флуктуаций центров масс всей системы и газовой и жидкой фаз в отдельности.

Основная трудность МД-моделирования связана, как известно, с тем, что для расчета сил, действующих на частицы, необходимо на каждом шаге по времени вычислять все расстояния между частицами. Число таких расстояний равно N(N - 1) /2, поэтому объем вычислений растет пропорционально Л'2. Наиболее удачный способ из- ; бежать квадратичного роста времени счета с увеличением числа ча- ••" стиц был предложен Верле [15]. Он состоит в определении матрицы ближайших соседей, которые только и учитываются при вычислении межчастичных расстояний. С этой целью вводятся-радиус обрезания взаимодействия гс и радиус горизонта гд . Частицы, находящиеся от данной ва рас£тчяции,1це превышающем гз ,*счйтаются ее ближайшими соседями. Во взаимодействии участвуют только частицы из сферы взаимодействия г;;- < ге. Частицы, расположенные между сферами с радиусами гс и г3 , образуют буфер, между которым и сферой взаимодействия происходит обмен частицами. Матрица «¡седей обновляется через каждые к шагов по времени. Чтобы даже бьгстрые частицы не могли преодолеть буфер за время кк , где к - шаг по времени, в

[15] было предложено выбрать радиус горизонта равным

r, = rc+sfc<u)/i, (7)

где (v) - среднеквадратичная скорость частиц и параметр s и 5,3 . В описанный классический алгоритм в данной работе были введены два новых элемента. Во-первых, частицы j , выходящие в течение времени kh за радиус горизонта »-й частицы, исключались из числа ее ближайших соседей до следующего обновления матрицы соседей. Эта процедура позволила ускорить работу программы примерно на 10% и дала возможность использовать формулу (8). Во-вторых, мы полагали, что г, не является постоянной величиной, как это считалось в (15), а зависит от скорости частицы и структуры ее ближайшего окружения. Мы полагали, что г,(i) увеличивается с ростом скорости частицы u¡ , уменьшается с увеличением числа частиц N¡ внутри радиуса обрезания i -й частицы и уменьшается по мере приближения к моменту обновления матрицы соседей. Указанным условиям удовлетворяет эмпирически подобранная зависимость

\ /•- rf (i, т) = Гс + (fc - т)(*)Ь (З + , + (8)

где (к — т) - число шагов по времени до обновления матрицы соседей, и,- - модуль скорости i -й частицы, С - эмпирическая постоянная. Путем предварительных расчетов подбиралось максимальное значение постоянной С , при котором не наблюдались случаи пробивания буфера внешними частицами. В наших расчетах было принято значение С = 0,88. Если использовать параметры нишей модели N(i) « 100, г> == 3,2 , то при т ~к/2, = (v) получается значение параметра! в « 2,2 . Это приводит к уменьшению числа соседей (при к — 16; (v) = 1/4 , /i = 1/33) no сравнению с (7) примерно в 1,5 раза. Эта оцецка позволяет заключать, что определение радиуса шризон-та по формуле (8) существенно повышает эффективность алгоритма по сравнению с исходным определением [15]. Более того, вследствие своей чувствительности к локальному значению нлотоости числа частиц формула (8) дает определенные преимущесгва при моделировании двухфазных систем. Отметим, что формула (8) не является оптимальной с точки зрения минимизации числа ближайших соседей. Этот вопрос требует дополнительного специального исследования.

5 В третьей главе приводятся результата исследования методом МЦ двухфазной ленпард-даонсовскоИ системы в состоянии термодинамического равновесия я при сташмиарном процесс« кспарспня о вакуум.

Система приводилась н равновесное состояние с заданной температурой, затем ее физические характеристики вычислялись путем усреднения вдоль фазовой траектории на поверхности постоянной энергии, т.е. моделировался микроканонический ансамбль. Как уже отмечалось, мы будем пользоваться МД-единицми. Во всех случаях число частиц в основной ячейке равнялось N = 8000 , ячейка имела размеры Ь х Ь х 2Ь с Ь = 25,1984. Шаг интегрирования составлял Н = 0,03125 , радиус обрезания гс = 3,2 . Укажем для справок размерные значения МД-единиц в случае аргона (атомная масса 40): единица скорости ьа = 1093,3 м/с, единица плотности числа частиц = 1,6825 г/см3. Единицы длины, времени и энергии были приведены выше, Некоторые результаты расчетов суммированы 6 Таблице 1.

Таблица 1. Параметры закрытых двухфазных систем

г 0.7250 ± 0.0058 0.7952 ± 0.0063

ni 0.8060 0.7722

п, 0.004886 0.009655

Hi -5.423 ± 0.5796 -5.146 ±0.6070

и, -0.05292 ±0.1364 -0.09212 ±0.1741

1.087 ±0.8874 1.190 ± 0.9721

ef 1.107 ±0.8991 1.232 ± 1.002

Т 0.897 ±0.007 0.991 ±0.008

ГЦ 0.7164 0.6582

п, 0.02367 0.04687 .

Ut -4.713 ±0.6463 -4.285 ± 0.6845

U, -0.2282 ¿0.2851 -0.4263 ±0.4013

(1 1.346 ± 1.099 1.486 ±1.214

ч 1.329 ±1.077 1 1.474 ±1.200

Для жидкой и газовой фаз (индексы lag соответственно) приведены значения плотности п/ и пд, средней потенциальной энергии на частицу щ и щ, средней кинетической энергии (/ а е/ при различных температурах (после знака ± даны значения средней квадратичной флуктуации). Как легко видеть, для флуктуаций кинетической энергии с высокой точностью выполняется равенство (((с) — е2)) = |(f2), что свидетельствует о термодинамическом равновесии в системе. Можно также отметить малую величину флуктуаций температуры ( « 0,8 %) вследствие большого числа частиц. Это оправдывает принятую постановку задачи, основанную на использовании

Рис. 1: Пространственные профили нлишичи числа ча< шц в шикнутиЙ cm кие при разных тешн-ратурах. Кршан 1 V = 0.725 ,2. Т - 0.795 ,3. Г = 0.897, Т- 0.991

NVE-ансамбля без искусственных П]>1»цолУ1> поддержания постоянной температуры.

Чтобы определить зависимость среднего знамения некоторой физической величины / от координаты г , МД-Яче11ка ризбштласьна 512 слоев толщиной ¿¡г = 2Л/512 и суммировались знамения величины / (которая зависит от координат н 'скоростей частиц) По всем частицам, попадающим и данный слой. Результат затем делился на полное число частиц, попавших и слой. Рассчитанные таким способом профили плотности числа частиц при разных температурах показаны на рис.1.

Как было описано н предыдущем разделе, вблизи дна ячейки на систему действует сила с hoiеццп алом (2), При Т = 0,725 постоянная полагались рннной 0,5 , при более высок и л температурах ft = 1 Действие зтого потенциала влияет на профиль платности только вблизи дна ячейки и не скашваетси на обымпых свойствах фаз

Рнс. 2: Пространственные профили средней потенциальной энергии иа одну частицу П и ее квадратичной флуктуации 6и в замкнутой системе при разных температурах. Обозначения кривых те же, что па рис.1.

Заметим, что в переходном слое между жидкостью и газом профиль плотности оказывается монотонным. Таким образом, наш расчет не подтверждает высказававшуюся ранее (см. [16]) гипотезу о существовании осцилляций плотности в переходном слое.

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

1'ис, 3: Средняя потенциальна* энергии на одну чьстицу а и ее флуктуации Ли кик функции плотности при речных температурах. Обозначения кривых le же, чю на рис.1.

На рис.3 средняя потенциальная энергия н ее флуктуация пршидс-11ы как функции плотности числи частиц. Криныс, cooriieii'i u.vioiiuic разным температурам, miuio ричлнчаютоя i» области идопиктеИ, tu-ощетствующей переходному слою. Поэтому для грубой оценки можно считать, что потенциальная чнерпш в переходном слое огм, функция локальной плотности.

Более детальную информацию «б энергии снят томов » жидкости, газе и переходном слое можно получить, рнсечнтыная но данным МД-моделирования функцию распределения потенциальной энергии частиц. Пример так!1« расчета предгтшюи на рис. 4 (Т » Ü.99J). Энергия частиц усреднялась по слоям толщиной 0.2 , параллельным плоскости ху. Координаты ¿ каждого из <v'|oea указаны н тибчице;

Рис. 4: Нормированные распределения потенциальной энергии р(и) для равновесной системы с Т = 0.991 , рассчитанный для различных слоев при переходе от жидкости (кривая 1) к газу (кривая 8). Описание п тексте.

1. - слой -15.0 < г < -14.0 5. - слой -4.0 < г < -3.8

2. - слой -6.5 < г < -6.3 6. - слой -3.5 < г < -3.2

3.-слой -5.5 < г <-5.3 7. - слой -3.0 < г < -2.8

4. - слой -4.5 < г < -4.3 8. - слой 3.0 < г < 18.0

Кривая 1 относится к однородной жидкости, кривая 8 - к газу, остальные кривые - к переходному слою. Для всех сечений функция распределения потенциальной энергии отлична от нуля при отрицательных энергиях. В газовой фазе она обнаруживает характерные особенности, связанные с образованием двухатомных и более крупных кластеров. Аналогичные особенности имеют распределения потенциальной энергии в сечеииях, близких к "газовому" краю переходного слоя (см. кривые 5 - 7 на рис. 4). Их можно интерпретировать как указание на зарождение кластеров внутри переходного слоя. Средняя потенциаль-

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

Важным источником информации о структуре переходной! слоя служит двухчастичная функция распределения п(г,-, Гу), Очевидно, что в однородной фазе она зависит только от межчастичного расстояния г,у = |г, — Гу|. В нашем случае, когда система неоднородна в направлении оси z, n(r,, Гу) является функцией трех переменных, в качестве которых удобно выбрать л -координаты двух частиц г;, гу и радиальное расстояние pfj — + у)2 между их проекци-

ями на плоскость ху. Двухчастичные функции распределения в указанных переменных были рассчитаны с использованием результант МД-моделирования. 4fo6u составить представление о том, как меняется функция n(p;y, г,-, Zj) при переходе от жидкой фазы к газовой, переходная область делилась илоскостямн, параллельными плоскости ху, на ряд слоев г, толщиной 0,2 и координата г; выбиралась внутри одного из этих слоев. Были исследованы распределения двух типов: 1 - радиальные, зависящие от переменной pij и рассчитанные нри условии |г, — гу| < 0.2 и 2 - аксиальные, зависящие от переменной (2j - Zi) и рассчитанные при условии /),у < 0,2 . Рас пределен и я обоих тшюн были построены также для объемных фаз. Расчеты производи лись при различных температурах or Т = 0,725 до Т = 0,991 . Ни рисунке 5 показаны радиальные распределения в различных слоях г,, координаты которых указаны в таблице:

1.- слой -15.0 < г < -чМ.О 5. - слой -4.5 < г <-4.3

2.- слой —6.5 < г < —Ç.3 С, - слой -4.0 < z < -3.8

3. - слой -6.0 < г < -3.8 7. - слой -3.5 < z < -3.2

4. - слой -5.5 < г < -5:3 8. слой 3.0 < z < 18.0

Видно, что при переходе от жидкой фазы К газоной в первую очередь разрушаются дальние координационные сферы. При ïiom положенно максимумов практически не меняется. Это означает, что наиболее нс-роятное расстояние между некоторой выделенной частицей и частицами » -ой координационной сферы стремится оПпгься неизменным в равновесной двухфазной системе. Таким образом, мри периоде m

Рис. 5: Парные радиальные распределения в равновесной системе при Т — 0.991 « зависимости от положения слоя г, при переходе от жидкости (кривая 1) к газу (кривая 8).

жидкости к газу проявляется' тенденция сохранения ближнего порядка •

Пример аксиального распределения показан на рис.б. Координаты слоев указаны в габлице:

1. - слой -15.0 < 2 < -14.0 4. - слой -4.0 < г < -3.8

2.-слой —6.5'< г < —6.3 5.-слой 3.0 < г < 18.0

3.-слой -5.5 < г <-5.3

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

Рис, 6: Парные аксиальные распределения в равновесной системе при < 0.2 в зависимости от положения слоя г", ( г = щ — 2,-) при Т = 0.991 . Квадраты на кривой 1 соответствуют радиальному распределевию в жидкости.

квадратами показаны радиальные распределения для жидкой фазы. На основе изучения парных корреляционных функций в переходной области между жидкостью и паром можно сделать следующее заключение о структуре поверхностного слоя жидкости. Вследствие сохранения ближнего порядка, локальная плотность (определенпая на масштабе 2-3 единиц длины) в межфазной области принимает только два значения, соответствующие плотностям объемных фаз - жидкости или газа. Это означает, что существует хорошо определенпая граница раздела между жидкостью и газом, толщина которой равна примерно 1-2 МД-единицам. Амплитуда флуктуацяй положения этой границы определяет средпий профиль плотности по оси г и среднюю толщину межфазпой области, которая увеличивается с ростом температуры и достигает примерно 8 МД-едишщ при Т « 0,9,

? = 0.2493 , 2. - д - 0.3526 , 3. - 9 = 0.7480

Мы сделали попытку применить для исследования флуктуации плотности в переходном слое формализм временных корреляционных функций. С этой целью в МД-ячейке были выделены три области, соответствующие объемным жидкой и газовой фазам и переходному слою. Для каждой из областей был проведен анализ флук-туаций плотности с волновыми векторами q/lm = ^(/,т,0) (I = 0,1,2,3;т = 0,1,2,3), лежащими в плоскости х,у и вычислен динамический структурный фактор. В момент времени * для каждой области вычислялся пространственный фурье-образ оператора плотности:

1 № ' п(я>*)= £ ехР('Чг>(0)> ,.

где ¿У, • число частиц замеченных в области в момент времени < , а

индекс j нумерует эти частицы. Временной ряд п^, <) записывался на постоянное запоминающее устройство. После окончания расчета с помощью быстрого преобразования Фурье проводилась оценка функции спектральной плотности стационарного процесса <), которая является искомым динамическим структурным фактором:

где = 0)п(я, I)) - автокорреляционная функция флук-

туаций плотности. Используя конечное преобразование Фурье это выражение можно переписать в виде

Здесь N - среднее число атомов в области. Вследствие статистической несостоятельности оценки спектральной плотности по единичной реализации случайного процесса [17] паи пришлось разбить полный временной интервал на отрезков с тем чтобы уменьшить ошибку оценки, пропорциональную 1 ¡-/Я}. При этом нами не использовались спектральные окна. Детали применяемой методики можно найти в [18]. Пример рассчитанного по описанной методике динамического структурного фактора показан на рис.7. При наибольших длинах волн и низких температурах (0.725,0.897) наблюдаются коллективные воз-бужения, скорость которых оказывается близкой к скорости звука в аргоне «1 ( 1000 м/с ) вблизи тройной точки.

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

где вектор q лежит в плоскости х, у. Развитый в [19] формализм допускает обобщение на динамический случай. Можно показать, что при

(10)

(И)

зависящих от времени флуктуациях плотности фурье-компонента высоты разделяющей поверхности Гиббса связана с фурье-компонентой флуктуаций плотности соотношением, подобным (11):

*о(ч, 0 -¡4гп((\, г, *)/(«/ — (12)

Если теперь подставить (12) в (10), то получим:

S(q,u) =(щ -щ) lim

1

г- ео 2nNT

J^ dt zq(<i, t) exp(—iwf)|

Таким образом, рассчитанный нами для переходного слоя динамический структурный фактор совпадает со спектром мощности флуктуаций разделяющей поверхности Гиббса. Для волновых векторов, лежащих в интервале 0,2493 < |к| < 0,899 , этот спектр представляет собой мощный низкочастотный шум. Исследование динамики разделяющей поверхности описанным выше методом требует расширения интервала волновых векторов.

Чтобы моделировать испарение жидкости в вакуум, потенциал на верхней грани МД-ячейки полагается равным нулю ("открытая крышка"). В этом случае атомы могут свободно покидать ячейку, и равновесные состояния не достигаются. Можно, однако, привести систему в стационарное неравновесное состояние с заданной постоянной температурой Т, на краю жидкой фазы. Это достигается переносом покинувших МД-ячейку частиц в придонную область ячейки с новыми скоростями, выбранными гак, чтобы сохранить полную кинетическую энергию системы и оставить равными нулю х - и у -компоненты полного импульса. В таком стационарном состоянии можно усреднять характеристики системы по времени и сравнивать результаты с характеристиками равновесной системы с Т — Т,.

Максимальная плотность потока атомов, испаряющихся с поиерх-ности конденсированного тела, дается формулой Герца {20): jm ~ п,^к[,Т/2лгп , где п, - плотность пара, находящегося в равновесии с жидкостью нри температуре Т • Отношение реальной плотности потока испаренных атомов j к ее максимальному значению j„t называется коэффициентом испарения. В следующей таблице приведены результаты расчетов потоков jm,j и коэффициента испарения а ~ j/jma* при различных температурах поверхности жидкости.

Рис. 8: Плотность числа частиц при разных температурах В зависимости от % : 1. - Т = 0.725 (равновесная система), 1', - Т, ез 0.725 (испарение в вакуум), 2. -Т = 0.897 (равновесная система), 2'. - Т, « 0.9 (испарение в вакуум)

т. Зтах 3 а

0.725 0.00023956 0.0001908 0,79646

0.795 0.00049577 0.0003932 0.79311

0.897 0.00129090 0.0010280 0.79636

Обращает па себя внимание независимость а от температуры. Возможно, что значение а ~ 0.8 является универсальным для простых жидкостей.

При моделировании испарения в вакуум наблюдается значительное изменеиие пространственных профилей плотности. Во-первых, сред-пяя плотность пара падает в несколько раз по сравиеяию с равновесной. Во-вторых, плотность объемных газовой и жидкой фаз становится неоднородной. Рис,8 показывает, что плотность жидкости Возрастает по линейному закону с ростом г . Из табл.2 видно, что плот-

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

Таблица 2. Параметры открытых двухфазных систем

Т. с1п1/(1г сЛ^/г/г

0.725 0.795 0.897 +0.00058 +0.00150 +0.00440 -0.000030 -0.000057 -0.000210 -0.005 -0.010 -0.034 -0.0020 -0.0038 -0.0110

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

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

Молекулярно-динамическое моделирование испарения в вакуум

рЫ • - ■ ■

1

• • • / • / ' . 1 V 1 | \ •

* 1 1 / • 1

________••"'1^/ "мм 114 1 Г 9]21 3\ 'Г<1 1 1 II II 1 1 1 » 1 I » 1 ^

Рис, И: Распределении частиц по скоростям р(и.) в зависимости от слоя г, при Т, « 0.897 :: + - переходной слой г, 6 (-6, -2) ,1. - переходной слой я, £ (-2,0), 2. - газовы» слой г.е (8,10), 3. - газовый слой г, € (22,24)

позволяет определить важнейшую характеристику процесса испарения — функцию распределения испареиных частиц по скоростям ^близи границы фаз. Зга функция, служащая граничным условием для газокинетических расчетов [21] обычно считается максвеллоаской зу1Я частиц, удаляющихся от поверхности (и/> 0), и равной нулю для движущихся к поверхности (ьж < 0). Сколько-нибудь серьезный теоретический вывод такого граничного условия нам не известен. Результаты расчета функций распределения частиц по поперечной скорости уя в объемных фазах и в переходном слое для различных температур , Достаточно близки к максвелловским с нулевой средней скоростью к температурой, несколько снижающейся при переходе от жидкости Щ- газу. Пространственные профили массовых скоростей и температур цохнзаиы на рис.0,10. Распре,цоленнв частин но продольной компонен-

0,5

0.4

0.3

0.2

0.1

0.0 •

0 5 10 15 2 20

Рис. 12: Доля частиц с отрицательной скоростью ввутр»слоев г, по мере удаления от границы жидкости: о - Г, = 0.725 , .+ - Т, = 6.897

те скорости V, носят,более сложный характер. Как видно из'рис.11 , внутри жидкости распределение максвелловское. Рисунок 12 показывает, что в переходном и газовом слоях почти монотонно уменьшается число частиц с отрицательными скоростями;-Однэ^о, даже у не] -ч н^й границы МД-ячейки доля частиц с отрицательными с'кброргями составляет около 8,1% для Т, = 0,725 и около 10,7% Д#я Т,"^'0,897.

Можно предположить существование двух механизмов образования ■: частиц с отрицательными скоростями в испарявшемся потоке.

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

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

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

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

Из рисунка 10 видно, что внутри МД-ячейки "температуры" рас' пределений по и по V, заметно различаются. Строго говоря, их даже нельзя сравнивать, поскольку функция распределения по V, существенно отличается от максвелловской. "Поперечная температура", соответствующая близким к максвелловскому распределениям по компонентам скорости ьХ1 ур , достаточно плавно меняется в переходном слое и газовой фазе. "Продольная температура", соответствующая распределению по 'ь, и имеющая смысл ширины распределения, резко падает в переходном слое (на масштабе 5-8 МД-единиц), а затем плавно меняется в газовой фазе, оставаясь в 1,5-2 раза меньше "поперечной температуры". Одновременно с резким изменением "продольной температуры" на том же масштабе толщины переходного слоя возникает конечная г -компонента массовой скорости. Доля частиц с отрицательной ? -компонентой скорости убывает с удалением от поверхности жидкости и должна стать пренебрежимо малой на расстояниях, много больших толщины переходного слоя. Таким образом, в случае малой плотности пара использование "полумаксвелловско-го" распределения в качестве граничного условия при решении гаю-кинетических задач представляется разумным ириближецием.

Четвертая глава посвящена применению молекулярно- динамического подхода к аналитическому рассмотрению фазового перехода кристалл- жидкость. ^ '

Известно, что кривая нлаьления нормальных веществ с хорошей точностью описывается эмпирическим уравнением Симона [22]:

' .•• Р^АТ^ + в / (13)

Автор диссертации предлагает, кок ему кажется, наиболее простой и

еще неизвестный способ вывода (13).

Рассмотрим систему N тел, взаимодействующих посредством потенциальной функции (7(г,..., г) = U(q), где i пробегает от 1 до 3N . Тогда классическая динамика будет определяться уравнениями Гамильтона (14) с И = К(р) +• U(g), где К(р) = р7/2т :

dt ~ dt dp,-' W

Предположим, что U(q) - однородная функция степени s :

ff (а®) = a'U(qi). Заметим, что K{p) - однородная функция степени два:

K(aPi) = a'Kipi)

Допустим существование некоторого решения (14): g¿ = q,(t),p¡ = Pi{t), такого, что фазовая траектория соответствует макросостоянию, представляющему собой равновесную многофазную систему. Тогда, можно сопоставить этому макросостоянию термодинамические переменные P,V,T. Поскольку мы используем NVE -ансамбль, то давление и температура будут флуктуировать, но в пределе N,V оо описанная процеду ра не вызывает сомнений. Отметим, что пристеночный потенциал легко вводится простым добавлением к U членов вза-имодествпя той Же степени s, но мы не будем усложнять написание формул так как конечный результат не зависит от этих членов.

Произведем следующее масштабное преобразование канонических переменных (в случае явного учета пристеночного потенциала необходимо дополнить (15) преобразованием координат стенок):

q' = aq, р' = bp, t' ~ et (15)

где а,Ь,с - константы, и потребуем, чтобы q' = q'(t') = aq(ct),p' = p\t') = bp(ci) также были решением (14). Для этого необходимо,чтобы преобразование (15) сохраняло форму уравнений (14). Используя однородность функций К и U можно легко получить требуемые условия:

б2 = а', <? = а1". (15а) '

где а - произвольная положительная константа.

Обратим внимание на главную особенность преобразования (15) -оно сохраняет геометрические отношения на фазовой траектории, тем

ч

25

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

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

У = а3Ъ .

- т=*'т° - у=ЧтГ (1б)

Здесь постоянная Больцмана % = 1 и (...) означает усреднение по времени. Показатель степени в (16) был известен и раньше, но для его получения были использованы другие методы. Замечательная особенность предлагаемого метода состоит в том, что любая физическая величина Г = <)) может быть представлена вдоль кривой плавления аналогично (16). Например, коэффициент самодиффузии связан с автокорреляционной функцией скорости формулой:

ОО

В = | ^ (у(0)у(<))Л . Применяя преобразование (15) получим зависимость О от температуры плавления:

О = Д,(ВД),/3-' (17)

•Интересно отметить, чтр для гипотетических кулоновских (или гравитационных) систем (в — ~1(> —2)). коэффициент самодиффузии с ростом температуры плавления падает, а для систем .с более резким отталкивательным потенциал^! (л < -2) наоборот расгет.

Прямым следствием масштабного преобразования (15) и формул (16) является отсутствие критической точки на кривой равновесия фаз в системах с однородной; потенциальной функцией.

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

ОСНОВНЫЕ РЕЗУЛЬТАТЫ .

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

2. Показано, что при переходе от жидкости к газу проявляется тенденция сохранения ближнего порядка. Вследствие этого локальная плотность в межфазной области принимает фактически только два значения, соответствующие плотностям объемных фаз -жидкости или газа. Это означает, что существует Хорошо определенная граница раздела между жидкостью и газом, толщина которой равна примерно 1-2МД-единицам. Амплитуда флуктуаций положения этой границы определяет средний профиль плотности и среднюю толщину межфазной области, которая увеличивается с ростом температуры и достигает примерно 8 МД-единиц при Т~ 0,9 МД-едйниц.

3. Испарение в вакуум незначительно влияет на структуру переходного слоя жидкость-пар при сравнении с равновесным случаем. Коэффициент испарения á не зависит от температуры поверхности жидкости в исследованном диапазоне, и, возможно, значение а ~ 0.8 является универсальным для простых жидкостей.

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

Основные результаты диссертации изложены в следующих работах:

1. Anisimov S.I., Barsukov A.V., Zhakhovskií V.V. Phase transition kinetic in metals subjected to laser pulses // Intern. Colloquium on Gas Kinetics in the Presence of Phase Transitions, Proc. of Euromech Colloquium N.248, Minsk, Byelorussia, Sept. 1991

2. Anisimov S.I., Barsukov A.V., Zhakhovskii V.V. Melting and evaporation of metals heated by picosecond laser pulses: changes of reflectivity // 3 rd. International Symposium on Subsecond Therinophysics, Yraz, Austria, September, 1992

3. Анисимов С.И., Жаховский B.B. Об испарении жидкости // Письма в ЖЭТФ - 1993. - т.57 - N.2 - с.91-94

4. Жаховский В.В., Анисимов С.И., Численное моделирование испарения жидкости методом молекулярной динамики // ЖЭТФ -1997. - (в печати)

5. Жаховский В.В. Кривая плавления систем с однородной потенциальной функцией // ЖЭТФ - 1994. - т.105 - N.6 - с.1615-1620

Литература

[1] Я.И.Френкель, Кинетическая теория жидкостей, Л.: Наука, 1975, 592 с.

[2] О.Кнаке, И.Н.Сгранский, УФН, т.68, N 2, с.261, 1959

[3] Д.Хирс, Г.Паунд, Испарение и конденсация, М.: Металл., 1966, 196 с.

[4] С.И.Анисимов, ЖЭТФ, т,54, с.338, 1968

[5] M.N.Kogan, Aun.Rev.Fluid Mech., v.5, p.383, 1973

[6] С.И.Анисимов, А.Х.Рахматулина, ЖЭТФ, т.64, с.869, 1973

[7J В.И.Жук, Изв. АН СССР, МЖГ, N.2, с.97, 1976

[8] T.Ytrehus, In:Rarefied Gas Dynamics, Ed. by J.L.Potter, NY, AIAA, 1977, p.1197

[9] C.Cercignani, In:Rarefied Gas Dynamics, Ed. by P.Camparque, Parie, CEA, 1979, v.l,p.l41 4

[10] C.Cercignani, The Boltzmann Equation aud its Applications, NY, Springer, 1988

[11] С.И.Анисимов, В.В.Жаховский, Письма в ЖЭТФ, т.57, с.91, 1993

[12] Л.Д.Ландау, Е.М.Лифшиц, Гидродинамика, М.: Наука, 1986, с.342

[13] Э.Хайрер, С.Нёрсетт, Г.Ваннер, Решение обыкновенных дифференциальных уравнений. Нежесткие задачи, М.: Мир, 1990.

[14] Д.В.Хеермаи, Методы компьютерного эксперимента в теоретической физике, М.: Наука, 1990.

[15] L. Verlet, Phys. Rev., v,159, p.98,1967

[16] И.З. Фишер, Статистическая теория жидкостей. М., Физматгиз, 1961,280 с.

[17] А.Н.Ширяев, Вероятность. М.: Наука, 1989, с.640

[18] Дж.Бендат, А.Пирсол, Прикладной анализ случайных данных, М.: Мир, 1989, с.540

[19] D.G.Triezenberg, R.Zwanzig, Phys.Rev.Letts, 1972, v.28, р. 1183

[20] Л.Д.Ландау, Е.М.Лифшиц, Статистическая физика. Часть 1, М.: Наука, 1976.

[21] Физическая кинетика и процессы переноса при фазовых превращениях. Под ред. С.И-Аниснмова, Минск.: Наука и техника, 1980.

Жаховский Василий Викторович

МОДЕЛИРОВАНИЕ ФАЗОВЫХ ПЕРЕХОДОВ ПЕРВОГО РОДА МЕТОДОМ МОЛЕКУЛЯРНОЙ ДИНАМИКИ

[22] С.М. Стишов, УФН, т.114, с.3-40,1974

Автореферат

Подписано к печати 15.10 96 Печать офсетная Тираж 100 экз.

Уч,-изд.л.2,оо Заказ N 130

Формат 60 х 84/16 Усл.печ.л. 1.86 Бесплатно

АП "Шанс". 127412, Москва, ул. Ижорская, 13/19

н

Введение

1 Метод молекулярной динамики и моделирование системы жидкость-пар

2 Метод МД для двухфазных систем

2.1 Интегрирование уравнений движения.

2.2 Математическая модель двухфазной системы

2.3 Вычислительный алгоритм.

3 МД моделирование системы жидкость-пар

3.1 Профили физических величин

3.2 Двухчастичные функции распределения

3.3 Распределение частиц по энергии.

3.4 Динамический структурный фактор.

3.5 Испарение в вакуум.

4 Кривая плавления

4.1 Вывод уравнения Симона.

4.2 Сравнение с экспериментом.

Актуальность. Проблема испарения конденсированного вещества относится к числу классических. Она представляет принципиальный интерес и важна для многих приложений в научных исследованиях и технологии. Интенсивному изучению этой проблемы в последние годы способствовала разработка эффективных методов генерации высоких плотностей энергии, основанных на использовании лазерного излучения и мощных пучков частиц. Полученный таким образом обширный экспериментальный материал обычно интерпретируется на основе весьма упрощенных теоретических моделей (см., например, [28, 31, 32]).

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

Целью настоящей работы является:

1. Создание высокоэффективного вычислительного алгоритма и программы моделирования двухфазных систем методом молекулярной динамики.

2. Численное моделирование термодинамически равновесной системы жидкость-пар и процесса испарения простой жидкости в вакуум.

3. Исследование структуры переходного слоя между объемными газовой и жидкой фазами, вычисление флуктуаций по-тециальной энергии и парных корреляционный функции в объемных фазах и переходном слое. Определение функции распределения атомов по скоростям при испарении в вакуум.

4. Выяснение роли флуктуаций энергии связи атомов в поверхностном слое в процессе испарения и определение механизма испарения простых жидкостей.

5. Определение вида кривых фазового равновесия (в частности кривой плавления) и зависимости коэффициентов переноса от температуры фазового перехода в системах с однородной потенциальной функцией.

Научная новизна работы состоит в следующем:

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

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

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

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

Результаты исследований обобщены в виде следующих положений, выносимых на защиту:

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

2. Показано, что при переходе от жидкости к газу проявляется тенденция сохранения ближнего порядка. Вследствие этого локальная плотность в межфазной области принимает фактически только два значения, соответствующие плотностям объемных фаз - жидкости или газа. Это означает, что существует хорошо определенная граница раздела между жидкостью и газом, толщина которой равна примерно 1-2 МД-единицам. Амплитуда флуктуаций положения этой границы определяет средний профиль плотности и среднюю толщину межфазной области, которая увеличивается с ростом температуры и достигает примерно 8 МД-единиц при Т ~ 0,9 МД-единиц.

3. Испарение в вакуум незначительно влияет на структуру переходного слоя жидкость-пар при сравнении с равновесным случаем. Коэффициент испарения а не зависит от температуры поверхности жидкости в исследованном диапазоне, и, возможно, значение а ~ 0.8 является универсальным для простых жидкостей.

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

Апробация работы. Результаты исследований были представлены на International Colloquium on Gas Kinetics in the Presence of Phase Transitions (Minsk, 1991), 3-rd. International Symposium on Subsecond Thermophysics, (Yraz, Austria, 1992), нескольких семинарах Теоретического отдела ИВТ РАН. По материалам диссертационной работы опубликовано 5 печатных работ [26],[27],[36],[39],[40].

Структура и объем диссертации. Диссертация состоит из введения, четырех глав и приложения, содержит 121 страниц машинописного текста, 57 рисунков на 57 страницах, 9 таблиц, 52 наименований использованной литературы.

1. Hoover W.G. Nonequilibrium molecular dynamics: the first 25 years. // Physica A. - 1993.- v. 194 - p.450-461.

2. Wainwrigt Т., Alder B.J. Molecular dynamics computations for the hard sphere system. // Nuovo cim.- 1958.- v.9 Suppl.l-p.116-131.

3. Alder B. J., Wainwrigt Т.Е. Studies in molecular dynamics. I. General method. // J.Chem.Phys. 1959.- v.31 - N.2 - p.459-466.

4. Alder B.J., Wainwrigt Т.Е. Studies in molecular dynamics. II. Behavior of a small number of elastic spheres. // J.Chem.Phys. -I960.- v.33 N.5 - p.1439-1451.

5. Alder B.J., Wainwrigt Т.Е. Enhancement of diffusion by vortexlike motion of classical hard particles. // J.Chem.Soc.Japan. -1969.- v.26 p.267-269.

6. Dymond J.H., Alder B.J. Vortex motion in fluids. // Ber. Bun-senges. Phys. Chem. 1971.-v.75-N.3-4-p.394-396.

7. Verlet L. Time-dependet correlations in simple fluids. // Ber. Bun-senges.Phys.Chem. -1971.-v.75-N.3-4 -p.389-392.

8. Rahman L. Correlation in the motion of atoms in liquid argon. // Phys. Rev.- 1964.- V.136A N.2- p.405-411.

9. Verlet L. Computer experiments on classical fluids. I. Thermody-namical properties of Lennard-Jones molecules. // Phys. Rev.-1967.- v.159 N.l- p.98-103.

10. Optiz A.C.L. Molecular dynamics investigation of a free surface of liquid argon. // Phys. Letters A 1974.- v.47 - p.439-103.

11. Rao M., Levesque D. Surface structure of a liquid film. //J. Chem. Phys. 1976.- v.65 - N.8 - p.3233-3236

12. Chapela G.A.,Saville G.,Thompson S.M.,Rowlinson J.S. Computer simulation of a gas-liquid surface. Part 1. // J. Chem. Soc. Faraday Trans.- 1977.- v.28 p. 1133-1144.

13. Nijmeijer M.J.P.,Bakker A.F.,Bruin C.,Sikkenk J.H. A molecular dynamics simulation of the Lennard-Jones liquid-vapor interface. // J.Chem.Phys.- 1988.- v.89 p.3789-3792.

14. Nijmeijer M.J.P., Bruin C., Bakker A.F., van Leeuwen J.M.J. Wetting and drying of an inert wall by a fluid in a molecular-dynamics simulation. // Phys.Rev.A- 1990.- v.42 N.10-p.6052-6059.

15. Haye M.J.,Bruin C. Molecular dynamics study of the curvature correction to the surface tension. // J.Chem.Phys 1994.- v.100- N.l p.556-559.

16. Yasuoka K., Matsumoto M., Kataoka Y. Evaporation and condensation at a liquid surface.I.Argon. // J.Chem.Phys.- 1994.- v.101- p.7904-7911.

17. Matsumoto M., Yasuoka K., Kataoka Y. Evaporation and condensation at a liquid surface.II.Methanol. // J.Chem.Phys.- 1994.-v.101 p.7912-7917.

18. Triezenberg D.G., Zwanzig R. Fluctuation theory of surface tension // Phys. Rev .Letts1972- v.28 p. 1183-1185.

19. Ubbelohde A.R. The Molten State of Matter. John Wiley and Sons Ltd.,1978.

20. Rowlinson J.S. The statistical mechanics of systems with steep intermolecular potentials. // Molecular Physics 1964. - v.8 -p.107-115

21. Marvin Ross Generalized Lindemann melting law. //Phys. Rev. -1969. v.184 - p.233-242

22. Hansen J.P. Phase transition of the Lenard-Jones system. II. High-temperature limit. // Phys.Rev.A -1970. v.2 - p.221

23. Hoover W.G, Gray S.G., Johnson K.W. Thermodynamic properties of the fluid and solid phases for inverse power potentials. //J. Chem. Phys.- 1971. v.55 - p.1128-1136

24. R.K.Crawford, Rare Gas Solids, vol.2 (edited by M.L.Klein and J.A.Venables). Academic Press, N.Y., 1977, p.663

25. Anisimov S.I., Barsukov A.Y., Zhakhovskii V.V. Phase transition kinetic in metals subjected to laser pulses // Intern. Colloquium on Gas Kinetics in the Presence of Phase Transitions, Proc. of Euromech Colloquium N.248, Minsk, Byelorussia, Sept. 1991

26. Anisimov S.I., Barsukov A.V., Zhakhovskii V.V. Melting and evaporation of metals heated by picosecond laser pulses: changes of reflectivity // 3 rd. International Symposium on Subsecond Ther-mophysics, Yraz, Austria, September, 1992

27. Френкель Я.И. Кинетическая теория жидкостей. Д.: Наука, 1975, - 592 с.

28. Фишер И.З. Статистическая теория жидкостей. М.: Физмат-гиз, 1961, - 280 с.

29. Роулинсон Дж.,Уидом Б. Молекулярная теория капиллярности. М.: Мир, 1986, - 376 с.

30. Кнаке О., Странский И.Н. Механизм испарения. // УФЫ -1959 т.68 - N 2 - с.261-305

31. Хирс Д., Паунд Г. Испарение и конденсация. М.: Металлургия, 1966, - 196 с.

32. Лабунцов Д. А. Неравновесные эффекты при испарении и конденсации. Сб. Парожидкостные потоки. - Минск, ИТМО им.А.В.Лыкова, 1977, с.6-33.

33. Муратова Т.М., Лабунцов Д.А. Кинетический анализ процессов испарения и конденсации. // ТВТ 1969 - т.7 - N 5 -с.959-967

34. Макашев Н.К. Испарение, конденсация и гетерогенные химические реакции при малых значениях числа Кнудсена // Ученые записки ЦАГИ, 1974 - т.5 - N 3 - с.49

35. Анисимов С.И., Жаховский В.В. Об испарении жидкости // Письма в ЖЭТФ 1993. - т.57 - N.2 - с.91-94

36. Анисимов С.И.,Рахматулина А.Х. Динамика расширения пара при испарении в вакуум // ЖЭТФ 1973. - т.64 - N.3 - с.869

37. Физическая кинетика и процессы переноса при фазовых превращениях. Под ред. С.И.Анисимова Минск.: Наука и техника, 1980, - 208 с.

38. Жаховский В.В., Анисимов С.И. Численное моделирование испарения жидкости методом молекулярной динамики // ЖЭТФ 1997. - (в печати)

39. Жаховский В.В. Кривая плавления систем с однородной потенциальной функцией // ЖЭТФ 1994. - т. 105 - N.6 - с. 16151620

40. Коган М.Н. Динамика разреженного газа. Кинетическая теория. М.: Наука, 1967, - 440 с.

41. Марч Н., Тоси М. Движение атомов жидкости. — М.: Металлургия, 1980, 296 с.

42. Ландау Л.Д., Лифшиц Е.М. Статистическая физика. Часть 1. М.: Наука, 1976, - 584 с.

43. Хайрер Э., Нёрсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. М.: Мир, 1990, - 512 с.

44. Хеерман Д.В. Методы компьютерного эксперимента в теоретической физике.- М.: Наука, 1990, 176 с.

45. Каплан И.Г. Введение в теорию межмолекулярных взаимодействий. М.: Наука, 1982, - 312 с.

46. Стишов С.М. Термодинамика плавления простых веществ. // УФН 1974. - т. 114 - с.3-40.- 121

47. Храпак А.Г. К теории плавления. // Письма в ЖЭТФ 1990. - т.51 - с.403-405.

48. Динамические свойства твердых тел и жидкостей. Исследования методом рассеяния нейтронов.(под ред. С.Лавси и Т.Шпрингера) М.: Мир, 1980, - 496 с.

49. Бендат Дж., Пирсол А. Прикладной анализ случайных данных. М.: Мир, 1989, - 540 с.

50. Ширяев А.Н. Вероятность. М.: Наука, 1989, - 640 с.

51. Справочник по специальным функциям с формулами, графиками и таблицами (под ред. М.Абрамовича и И.Стиган). М.: Наука, 1979, - 830 с.