Явления тепло- и массопереноса в пылевой плазме

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

Введение

Пылевая плазма (dusty plasma) представляет собой ионизированный газ, содержащий заряженные частицы конденсированного вещества. Для обозначения таких систем используются также термины “комплексная плазма” (complex plasma), “коллоидная плазма” (colloidal plasma), “плазма с конденсированной дисперсной вазой” (КДФ).

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

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

В лабораторных условиях пылевая плазма была впервые обнаружена Лэнгмюром в 1920-х годах [10]. Не смотря на это, ее активное изучение началось лишь в последние десятилетия в связи с развитием приложений, таких, как электрофизика и электродинамика продуктов сгорания ракетных топлив, электрофизика рабочего тела магнитогидродинамических (МГД) генераторов на твердом топливе и прочие. В конце 80-х годов исследования сосредоточились в основном на изучении зарядки пыли, распространения электромагнитных волн, их затухания и неустойчивости, как правило, применительно к пылевой плазме в космосе. В начале 90-х годов рост интереса к пылевой плазме был связан с широким использованием технологий плазменного напыления и травления в микроэлектронике, а так же при производстве тонких пленок и наночастиц. В середине 90-х годов удалось наблюдать формирование кристаллических структур в различных типах плазмы. Все это послужило толчком к бурному росту исследований в данной области, который продолжается и по сей день.

Среди современных направлений исследований в области пылевой плазмы основными являются:

1) образование упорядоченных структур, кристаллизация и фазовые переходы в системе пылевых частиц в различных типах плазмы;

2) элементарные процессы в пылевой плазме: зарядка пыли; взаимодействие между частицами в плазме; явления тепло- и массопереноса; внешние силы, действующие на пылевые частицы;

3) линейные и нелинейные волны в пылевой плазме, их динамика, затухание и неустойчивости.

Данная работа относится ко второму направлению. Рассмотрим результаты исследования процесса диффузии, полученные ранее.

Впервые формула для коэффициента диффузии была выведена Эйнштейном в 1905 году. Она связывала средний квадрат смещения (СКС) частицы и времени линейным соотношением , где - размерность пространства. Позже было обнаружено, что во многих случаях данная зависимость не является линейной. Процессы, в которых не выполняется формула Эйнштейна, называются процессами аномальной диффузии. Аномальная диффузия наблюдается как в жидкостях и газах [2], так и в плазме [3]. В качестве простого примера рассмотрим идеальный газ. Его молекулы взаимодействуют между собой как упругие сферы (то есть только при столкновениях), таким образом, рассмотрев малые времена, на которых частицы еще не успели столкнуться с другими, мы получим, что их движение является равномерным ,то есть зависимость СКС от времени является квадратичной. Данный пример показывает, что аномальная диффузия будет наблюдаться на достаточно малых временах, однако на практике она может наблюдаться и на бомльших временах, что показано, в частности, в работе. В ряде статей [3-5] была рассмотрена диффузия (как классическая, так и аномальная) в пылевой плазме. В них предоставлены как экспериментальные, так и теоретические результаты, полученные путем моделирования. При моделировании пылевой плазмы в них не были учтены флуктуация зарядов частиц и/или особенности приэлектродного слоя, что позволяется надеяться на улучшение полученных в них результатов.

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

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

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

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

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



1. Основные силы, действующие в пылевой плазме

1) Сила тяжести.

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

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

2)Сила трения со стороны нейтрального газа.

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

,

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

3) Случайные соударения с молекулами нейтрального газа.

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

, (1)

где -- нормально распределенная случайная сила.

4) Взаимодействие между пылевыми частицами в плазме.

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

Как показано в [1] в случае в изотропной плазме для потенциала вокруг уединенной сферической частицы с зарядом можно использовать выражение

(2)

Однако так же указывается, что рассматриваемый равновесный случай редко реализуется в плазме. Это объясняется тем, что из-за поглощения электронов и ионов пылевыми частицами, для существования плазмы при наличии пыли необходимы постоянные источники ионизации, поддерживаемые подводом энергии в разряд. То есть система является открытой. При этом распределения электронов и ионов в окрестности частицы оказываются неравновесными. Вопрос о том, в каких случаях можно использовать потенциал (2) был рассмотрен, в частности, в работах [67, 88 из Фортова]. Основные результаты этих работ можно сформулировать следующим образом: для частиц малого размера () на расстояниях, вплоть до существенно превышающих потенциал (2) хорошо аппроксимирует результаты, полученные в результате численных расчетов.

5) Электростатическая сила.

При наличии в плазме электрического поля напряженности на заряженную частицу действует сила

.

Введя эффективную величину поля

.

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

Для пылевой плазмы характерны условия, при которых , таким образом,

.

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



, (3)

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

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

. (4)

2. Эффекты, наблюдаемые в пылевой плазме

Рассмотрим эффекты, вносящие существенный вклад в поведение пылевых частиц.

1) Зарядка пылевых частиц.

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

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

2) Флуктуации заряда пылевых частиц.

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

, (5)

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

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

, (6)

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

Более того, даже в изотропной пространственно-однородной невозмущенной плазме заряд частиц будет колебаться около своего равновесного значения. Это связано с тем, что электроны и ионы поглощаются пылевыми частицами не непрерывно, а дискретно, в случайные моменты времени и в случайной последовательности. Учитывая это, можно выразить зависимость заряда не только от положения частицы, но и от времени:

, (7)

где нормированные флуктуации заряда описываются корреляционной функцией [6]:

, (8)

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

,

где -- зарядовое число.

3) Прочие эффекты.

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



3. Характеристики системы пылевых частиц

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

, (9)

где - концентрация пылевых частиц;

2) также часто рассматривают нормированный параметр неидеальности , который связан с параметром соотношением

; (10)

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

; (11)

4) соответствующая ей нормированная плазменная частота

; (12)

5) с последними связан параметр, который характеризует отношение силы взаимодействия между частицами и силы трения и выражается как



. (13)

Рассмотрим поведение параметра неидеальности:

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

b) также, из предыдущего следует, что Г возрастает, при увеличении заряда частиц или радиуса дебаевской экранировки и уменьшается при увеличении температуры системы.

4. Диффузия в пылевой плазме

Есть два основных метода, по которым рассчитывают коэффициент диффузии:

1) Формула Эйнштейна:

, (14)

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

2) Формула Грина-Кубо:

. (15)

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

Такое различие получаемых результатов обусловлено тем, что при выводе формулы Эйнштейна предполагалось, что автокорреляционная функция скорости (АКФС) убывает экспоненциально с ростом времени. Это предположение не всегда верно, однако при численном вычислении интеграла в формуле Грина-Кубо верхний предел берется достаточно большим, чтобы “хвост” АКФС давал малый вклад в итоговое значение. То есть в таком случае отклонение характера затухания АКФС от экспоненциального уже не играет важной роли.

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



, (16)

, (17)

для диффузии в вертикальном направлении, и

, (18)

, (19)

для диффузии в горизонтальном направлении.

Явление аномальной диффузии рассматривается, в частности, в работе [5]. В ней рассмотрены влияния различных характеристик системы на тип наблюдаемой диффузии. Так сначала рассматривается случай слабой силы трения (): для фиксированного и различных , меняющихся от 1 (случай неидеального газа) до 180 (точка плавления), строятся зависимости .

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

Что бы исследовать влияние скорости затухания на тип диффузии авторы строят графики для четырех значений Г (10, 50, 100 и 150) и различных .

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

Наконец было рассмотрено влияние коэффициента на процесс диффузии.

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

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

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

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

В этой статье было исследовано, влияние характеристик системы на характер наблюдаемой диффузии, а их влияние на значение коэффициента диффузии рассмотрено не было. Это влияние в сильно неидеальных системах было исследовано в статье [4].

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

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

. (20)

Действительно: при меняющемся от 50 до 102 зависимость, приведенная на рисунке 5(b) практически линейная, более того, зависимости, отвечающие различным значениям и , убывают приблизительно с одной скоростью (то есть ).

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

Далее авторы уточняют данную зависимость, приводя, в итоге, две формулы для предельных случаев:

, при (21)

, при . (22)

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

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

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

,

где -- количество частиц.

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

Говоря об этой работе, стоит так же отметить два положительных отличия от рассмотренных ранее статей.

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

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

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

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

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

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

Как видно, СКС вдоль оси всегда превосходит смещение вдоль оси.

5. Метод молекулярной динамики

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

, (23)

где - взаимодействующая сил, действующих на -ую частицу, а - число частиц.

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

Впервые метод молекулярной динамики был применен в 1957 году, в работе Олдера и Вэйнрайта [7], в которой на ЭВМ были проинтегрированы классические уравнения движения твердых сфер. С тех пор было проведено множество работ и исследований, опирающихся на ММД, а результаты, полученные в них, согласовываются с экспериментальными данными. Однако также возник и ряд вопросов об обоснованности применения данного метода, вызванный следующими наблюдениями.

1) Полная энергия в ММД не сохраняется, а флуктуирует, сохраняясь только в среднем. Это противоречит закону сохранения энергии, который должен выполняться для уравнений классической механики. Более того, при применении некоторых вычислительных схем наблюдался даже дрейф полной энергии.

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

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

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

3) Было обнаружено, что ММД прекрасно описывает необратимые процессы в полном согласовании с экспериментальными данными. То есть в методе, основанном на решении обратимых уравнений Ньютона, возникает необратимость.

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

Продемонстрируем это на примере сложения.

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

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

Не смотря на это, результаты, получаемые ММД, хорошо согласуются с экспериментальными данными, это объясняется тем, что в ММД практический интерес представляет не отдельная траектория, а результат усреднения по полному распределению траекторий. А этот результат оказывается устойчивым [8].

Так же в этой статье авторами делается важный вывод:

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

, (24)

для схем с явным вычислением скорости на каждом шаге или



(25)

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

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

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

6. Скоростная схема Верле

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

1) из известных на начало шага данных (, , где -- номер частицы, а -- номер шага) рассчитываются скорости частиц на полушаге:

,

где -- шаг численного интегрирования;

2) используя известные на начало шага и вычисленные находятся радиус-векторы :



;

3) полученные ранее данные подставляются в исходную систему уравнений, откуда вычисляются ускорения частиц в конце шага:

;

4) используя полученные ускорения находятся скорости в конце шага:

.

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

1) Решение задачи Коши для обыкновенного дифференциального уравнения.

Рассмотрим случай движения тела, на которое действует только сила трения.

Решая данное уравнение относительно получим



,

что с учетом начальных условий дает нам результат

.

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

Погрешность мала и не возрастает больше определенного значения.

2) Параметрический резонанс.

В качестве более сложного теста рассмотрим уравнение

(26)

и, решая его численно, найдем первые две зоны параметрического резонанса.

Кратко опишем способ построения данных областей:

а) Задается фиксированное значение , в данном случае оно было равно .

b) Для начального значения численно решается уравнение (26) при различных значениях отношения . Точки, в которых наблюдается параметрический резонанс, отмечаются на графике. Так получается одна горизонтальная линия на приведенном графике.

c) Значение изменяется, и находится следующая линия из зоны параметрического резонанса.

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

При таком подходе могут быть ошибки с определением того, является ли случай резонансным или нет. Так за 100 секунд амплитуда может возрасти в 39 раз и продолжать расти потом, однако случай ошибочно будет определен как нерезонансный. Эту погрешность определения можно уменьшить, проверяя, возросла ли амплитуда не в 40, а, например, в 10 раз, но в таком случае возникнет другая проблема -- амплитуда может возрасти в 10 раз, а потом начать уменьшаться (такой процесс называется биениями), то есть случай будет ошибочно признан резонансным. Биения наблюдаются на границе резонансной зоны и для данной задачи отношение амплитуд равное 40 при биениях не достигается. Таким образом, при автоматическом определении резонанса могут возникать погрешности, не связанные с точностью вычислительной схемы, однако оба явления, описанных выше (биения и слабый резонанс), наблюдаются вблизи границы области резонанса. При исследовании данной задачи, было проверено, что расстояние от точек, в которых могут наблюдаться эти явления, до границы области резонанса заметно меньше, чем шаг, с которым мы изменяем . То есть ошибки, связанные с определением того, является ли случай резонансным, не вносят дополнительной погрешности.

Полученные в итоге результаты согласуются, например, с результатами, приведенными в [9].

3) Термостат Ланжевена.

Наконец рассмотрим уравнение



, (27)

где -- случайная сила, удовлетворяющая следующим условиям:

,

,

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

Это уравнение описывает термостат Ланжевена, часто использующийся в работах по моделированию физических процессов для приведения системы к нужной температуре. Физически, такое уравнение движения описывает движение частицы, на которую действуют две силы: сила трения и сила упругих соударений с молекулами окружающего нейтрального газа (см. §1 глава1). Таким образом, константа выбирается в зависимости от типа и плотности нейтрального газа, а константа выбирается так, что бы данное уравнение описывало соударения с молекулами нейтрального газа, находящегося при данной температуре.

Решим данное уравнение относительно скорости

.

Пусть в начальный момент времени скорость частицы равна . Решение данного уравнения имеет вид

,



откуда получаем

,

.

Тогда для скорости получаем следующее выражение

.

Из него следую два важных соотношения:

, (28)

. (29)

Таким образом, среднее значение скорости стремится к нулю с течением времени, а средний квадрат скорости со временем стремится к постоянному значению, равному

.

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

, (30)



где -- постоянная Больцмана.

Теперь вернемся к численному решению уравнения.

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

, (31)

где -- шаг численного интегрирования. Таким образом, мы будем решать уравнение

. (32)

График квадрата скорости, полученного при численном решении данного уравнения, и его теоретическое значение.

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

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

7. Построение модели пылевой плазмы

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

1) Сила тяжести.

(33)

2) Взаимодействие с нейтральным газом.

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

, (34)

Где случайная сила подчиняется стандартному нормальному закону распределения.

3) Сила взаимодействия частиц.

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

. (37)



То есть соответствующая сила, действующая на -ую частицу, равна

, (36)

.

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

, (38)

соответственно электростатическая сила будет равна

.

При этом начало координат выбирается так, что при электростатическая сила и сила тяжести уравновешиваются, то есть

. (39)

5) Потенциальная ловушка.

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

,

. (40)



Еще раз отметим, что данная сила действует только в горизонтальной плоскости, то есть под в данном случае подразумевается проекция радиус-вектора на плоскость .

6) Изменение заряда пылевой частицы из-за смещения по вертикали.

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

можно пренебречь членами степени большей двух. То есть

. (41)

7) Случайные флуктуации заряда пылевой частицы.

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

a) Генерируем число поглощений заряженных частиц , произошедших за рассматриваемый промежуток времени. Это число подчиняется пуассоновскому закону распределения с интенсивностью .

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



и будет равняться нормированному значению изменения заряда за данный отрезок времени.

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

, (42)

где подчиняется стандартному нормальному закону распределения.

Учитывая пункты 4, 6 и 7 электростатическую силу можно записать в виде

(43)

Таким образом, движение системы пылевых частиц описывается системой уравнений

. (44)

8. Программный комплекс

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

1) Генерация нормально распределенной случайной величины.

Генератор псевдослучайных чисел в языке С++ выдает последовательность чисел, которая строится по переданному ему значению (random seed). То есть если будет передано одно и то же значение, то и последовательности псевдослучайных чисел будут одинаковыми. Однако если в качестве random seed передавать каждый раз новое число, то и последовательности будут различными. Для этого хорошо подходит функция time, которая при вызове с параметром NULL возвращает время, прошедшее с 1-ого января 1970 года, в секундах.

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

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

b) Первый вариант преобразования Бокса-Мюллера.

Генерируются две независимые случайные величины и , распределенные равномерно на промежутке . Найденные по формулам

(45)

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

с) Второй вариант преобразования Бокса-Мюллера.

Генерируются две независимые величины и , равномерно распределенные на отрезке . Если значение лежит в промежутке , то по формулам

(46)

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

На рисунке 13 приведены плотности вероятности величин, полученных с использованием центральной предельной теоремой и с использованием преобразования Бокса-Мюллера (два его варианта дают неразличимые результаты, поэтому приводится график только для второго варианта)

График расходится на две линии, это означает несимметричность функции распределения относительно нуля. На рис. 13(b) это расхождение гораздо меньше, особенно учитывая то, что он построен для гораздо более широкого диапазона значений случайной величины. Таким образом, хоть оба распределения и близки к нормальному, преобразования Бокса-Мюллера дают лучший результат. Более того, очевидно, что первый способ проигрывает им и в скорости.

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

2) Усреднение результатов.

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

a) после приведения системы к равновесию, рассматривается один отрезок по времени, превосходящий по длине тот, на котором мы хотим получить результат;

b) смещения находятся на пересекающихся отрезках интересующей нас длины:

c) усреднение производится по полученным .

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

Ускорение работы программы, вызываемое переходом к такому способу усреднения, настолько очевидно, что оцениваться не будет. В качестве проверки возможности такого перехода было проведено сравнение коэффициентов диффузии, посчитанных обоими способами, для частицы, движение которой описывается термостатом Ланжевена (время динамической памяти для такой системы, как видно из рисунка 8, около одной секунды). Результаты усреднения двумя способами мало различаются, эти различия не имеют выраженной зависимости от температуры и связаны с тем, что даже усредненные значения колеблются от запуска к запуску. Так же видно, что коэффициент диффузии возрастает при увеличении температуры.

9. Проверка возможности использовать метод молекулярной динамики

1) Так как для описания системы пылевых частиц используются уравнения классической механики, надо убедиться, что в системе не возникает квантовых эффектов. Это требование, безусловно, выполнено: мы рассматриваем движение только пылевых частиц (движение электронов и ионов нас интересует только в смысле их потоков на поверхность частицы), радиус которых равен . Конечно, такое рассмотрение является приближенным, однако на то, что бы описать еще и движение электронов и ионов (не важно, уравнениями классической или квантовой физики) нужны вычислительные мощности, недоступные на данный момент.

2) Проверим, что отсутствует дрейф полной энергии, то есть возможность использовать скоростную схему Верле.

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

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

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



Заключение

В рамках данной работы были достигнуты следующие результаты:

1) Изучен метод молекулярной динамики, условия его применимости, найдено время динамической памяти исследуемой системы.

2) Исследован термостат Ланжевена, установлен способ генерации случайной силы, при котором термостат приводит систему к нужной температуре. Для различных способов усреднения результатов (по независимым отрезкам и по времени) получены коэффициенты диффузии при различных температурах. Данные подтверждают рост коэффициента диффузии при увеличении температуры и равносильность рассмотренных способов усреднения.

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

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

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

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

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

Размещено на Allbest.ru

...