«Д. А. Медведев, А. Л. Куперштох, Э. Р. Прууэл, Н. П. Сатонкина, Д. И. Карпов МОДЕЛИРОВАНИЕ ФИЗИЧЕСКИХ ПРОЦЕССОВ И ЯВЛЕНИЙ НА ПК Учебное пособие НОВОСИБИРСК 2010 УДК 53.072 ББК 32.973 М 42 Медведев Д. А., Куперштох А. ...»

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РФ

НОВОСИБИРСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

СПЕЦИАЛИЗИРОВАННЫЙ УЧЕБНО-НАУЧНЫЙ ЦЕНТР НГУ

Д. А. Медведев, А. Л. Куперштох, Э. Р. Прууэл,

Н. П. Сатонкина, Д. И. Карпов

МОДЕЛИРОВАНИЕ ФИЗИЧЕСКИХ

ПРОЦЕССОВ И ЯВЛЕНИЙ НА ПК

Учебное пособие

НОВОСИБИРСК

2010 УДК 53.072 ББК 32.973 М 42 Медведев Д. А., Куперштох А. Л., Прууэл Э. Р., Сатонкина Н. П., Карпов Д. И. Моделирование физических процессов и явлений на ПК: Учеб. пособие / Новосибирск: Новосиб. гос. ун-т., 2010.

101 с.

ISBN 978-5-94356-933-3 Учебное пособие написано на основе преподавания одноименного спецкурса в Физико-математической школе при Новосибирском государственном университете (ныне Специализированный учебно-научный центр НГУ). Основное внимание уделено физическому содержанию задач, но детально разобран также ряд вычислительных алгоритмов, необходимых для моделирования.

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

Рецензент докт. физ.-мат. наук, профессор А. П. Ершов ISBN 978-5-94356-933- c Новосибирский государственный университет, c СУНЦ НГУ, c Медведев Д. А., Куперштох А. Л., Прууэл Э. Р., Сатонкина Н. П., Карпов Д. И., Оглавление Введение 1. Обязательные задания 1.1. Сложение колебаний..................... Задания............................. 1.2. Фигуры Лиссажу....................... Задания............................. 1.3. Движение маятника...................... 1.3.1. Безразмерные переменные............... 1.3.2. Фазовые диаграммы.................. 1.3.3. * Физический маятник.................. Задания............................. 1.4. Движение в центральном поле (планета).......... Задания............................. 1.5. Случайные блуждания.................... 1.5.1. Одномерные случайные блуждания......... 1.5.2. Средние значения................... 1.5.3. Функция распределения................ 1.5.4. Многомерные случайные блуждания......... 1.5.5. Случайные блуждания во внешнем поле.

Изотермическая атмосфера.............. Задания............................. 2. Задания для самостоятельной работы 2.1. Неравновесная агрегация, фракталы............ 2.1.1. Агрегация, ограниченная диффузией........ 2.1.2. Фрактальная размерность............... 2.1.3. Бессеточная модель.................. 2.1.4. Химически-ограниченная агрегация......... 2.1.5. Баллистическая агрегация............... 2.1.6. Кластер–кластерная агрегация............ 4 Оглавление Задания............................. 2.2. Электрический пробой.................... 2.2.1. Вычисление электрического потенциала....... 2.2.2. Критерии и модели роста стримеров......... Задания............................. 2.3. Молекулярная динамика................... 2.4.4. * Модель LBE с внешними силами и фазовыми 2.6. Теплопроводность, детерминированное горение...... 2.8.2. Ангармоническая цепочка.

Введение В течение многих лет на физическом факультете НГУ проводится компьютерный практикум ЧМФП (численное моделирование физических процессов) [1–3]. Практикум был создан по инициативе Л. М. Альтшуля, основной вклад в его развитие был сделан Г. Л. Коткиным. В СУНЦ НГУ этот практикум был существенно переработан, а также дополнен целым рядом новых постановок.

Это пособие представляет собой описание спецкурса Моделирование физических процессов и явлений на ПК для СУНЦ НГУ. Мы постарались уделить основное внимание физическому содержанию задач. Предполагается, что учащиеся знакомы с языком программирования Pascal, на котором выполнены программы.

Для того, чтобы основное внимание сосредоточить на собственно численном эксперименте, а не на программировании ввода-вывода и представления результатов, в курсе используется пакет графических программ и процедур обработки данных MPP (Modeling of Physical Processes), разработанный В. Г. Тупицыным, Т. Э. Овчинниковой и А. Л. Куперштохом. При этом на экране компьютера имеется графическое окно, в котором отображается поведение моделируемой системы, и таблица, где можно наблюдать и интерактивно изменять значения параметров и переменных. Описание графического интерфейса и основных процедур содержится в файле instrmp.txt. Некоторые дополнительные программы описаны в файле instr2.txt.

Спецкурс состоит из двух частей. На первом этапе учащиеся выполняют несколько обязательных заданий (обычно 4). Задания рекомендуется выполнять по порядку. К каждому заданию есть программа-заготовка.

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

На втором этапе предполагается самостоятельная работа над выбранным проектом. Некоторые темы для исследования описаны в главе 2.

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

В конце пособия приведен список рекомендуемой литературы. Во многом изложение следует книге [4]. Интересующимся алгоритмами расчетов рекомендуем книгу [5] замечательное пособие по методам вычислений. Остальные работы могут быть полезны желающим более подробно исследовать ту или иную проблему.

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

Обновленный вариант этой книги доступен в электронном виде по адресу http://ancient.hydro.nsc.ru/MPP_Specourse/mpp.pdf.

Глава 1.

Обязательные задания 1.1. Сложение колебаний Программа beats.pas.

Многие процессы в природе периодически повторяются. Математически это можно выразить в виде зависимости от времени f (t + T ) = f (t).

Величину T называют периодом функции. Часто используют также величину = 1/T, называемую частотой, и = 2 = 2/T круговую частоту.

Простейшая нетривиальная периодическая функция f (t) = = A sin(t + 0 ) описывает гармонические колебания. Величина A амплитуда колебаний, 0 начальная фаза.

В данном задании рассматривается сложение различных гармонических колебаний. Рассмотрим вначале сумму двух колебаний y(t) = = A1 sin(1 t)+A2 sin(2 t). Такая сумма может возникать, например, при сложении двух волн (волн на поверхности воды, звуковых, радиоволн) в случае, когда амплитуды их не слишком велики, так что выполняется принцип суперпозиции.

Результаты расчетов удобно представлять в графическом виде. Построим график функции y(t) = A1 sin(1 t) (чтобы ее получить, достаточно задать в таблице A2 = 0). Все знают, что график это линия, состоящая из бесконечного числа точек значений функции при непрерывном изменении аргумента. Понятно, что на практике мы можем вычислить значение функции только в конечном числе точек, например, в точках t = nt, n = 1, 2,... Как выбрать шаг t? Хотелось бы выбрать его побольше, чтобы за данное число операций просмотреть график на как можно большем интервале времени. Очевидно, однако, что слишком сильно увеличивать t тоже нельзя (начиная с t T полученная картина будет иметь мало общего с графиком функции y(t)).

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

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

При сложении двух колебаний с близкими частотами наблюдаются биения на экране видны высокочастотные колебания, амплитуда которых периодически меняется, но с меньшей частотой. Наглядно объяснить, что происходит, можно с помощью векторной диаграммы. Величину A sin(t) можно представить как проекцию на ось Y стрелки длиной A, вращающейся вокруг начала координат с угловой скоростью. Тогда сумма двух колебаний с частотами 1 и 2 это сумма проекций двух стрелок, или проекция суммы двух векторов. Если амплитуды одинаковы, суммарный вектор направлен под углом (1 + 2 )t/2 (вращается с угловой скоростью (1 + 2 )/2, дает высокочастотное заполнение). Длина его меняется в соответствии с изменением угла между стрелками со скоростью 1 2 и отвечает за амплитуду огибающей.

Можно посмотреть, что получится при сложении большого числа гармонических колебаний, если их частоты кратны основной k = k. Такая сумма называется рядом Фурье, ей можно аппроксимировать любую периодическую функцию. Слагаемые с частотами, кратными основной частоте, называются гармониками, например, слагаемое с частотой 9 девятая гармоника.

Задания 1. Выключите одно из колебаний, а частоту второго установите = 21, не меняя при этом шага по времени (t = 0,15). Объясните получившийся график (две пересекающиеся синусоиды). Можно попробовать другие частоты. Что следует предпринять, чтобы высокочастотные колебания отображались правильно?

2. Верните параметры в начальное состояние (можно просто перезапустить программу). На графике видны высокочастотные колебания, амплитуда которых изменяется так, что наблюдается низкочастотная огибающая. Каково отношение их частот? Как эти частоты получаются из 1 = 1 и 2 = 1,2?

3. Измените программу, чтобы складывались три колебания, частоты и амплитуды которых можно было бы задавать в таблице. Какой график получается при частотах 1 = 1, 2 = 3, 3 = 5 и амплитудах A1 = 1, A2 = 1/3, A3 = 1/5?

4. Получите график суммы N колебаний, число которых можно было бы задавать интерактивно. Интересны следующие зависимости амплитуды от частоты:

В двух первых случаях на разрывах функции видны вертикальные выбросы, высота которых не уменьшается при увеличении числа слагаемых и составляет около 18 % от высоты ступеньки (явление Гиббса, о нем можно прочитать в [6]).

1.2. Фигуры Лиссажу Программа liss.pas.

Подадим на вход Y осциллографа один гармонический сигнал y(t) = = A1 cos(1 t + 1 ), на вход X другой x(t) = A2 cos(2 t + 2 ). Можно убрать фазу в выражении для x, если сдвинуть начало отсчета времени на t = 2 /2, получим x(t) = A2 cos(2 t), y(t) = A1 cos(1 t + ), где = 1 1 2 /2. На экране мы увидим кривую линию, которая задана в параметрическом виде, то есть x и y зависят от параметра t.

Если отношение частот 1 и 2 рациональное число, то эта кривая будет периодически повторяться, и мы увидим неподвижную картинку фигуру Лиссажу. Иначе кривая линия в конце концов заполнит весь прямоугольник |y| A1, |x| A2. При некоторых значениях отношения частот и амплитуд получаются простые геометрические фигуры (прямые, эллипсы и т. д.). В этом задании Вам предлагается исследовать вид фигур Лиссажу при различных параметрах сигналов.

Задания 1. Не меняя ничего в программе (изменяя только значения амплитуд и частот, выведенные в таблицу), получите на экране горизонтальную прямую, вертикальный отрезок, отрезки под углом ±45.

2. Измените программу так, чтобы получить круг, восьмерку, знак бесконечности.

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

4. Задайте частоту по одной из координат равной, а по другой +, где (например, = /1000). Что происходит с графиком? В этом случае для вывода удобно использовать процедуру Orbit, пример ее применения есть в программе planet.pas.

1.3. Движение маятника Программа pendul.pas.

Математический маятник представляет собой точечное тело массы m, подвешенное на невесомом нерастяжимом стержне длиной l в поле тяжести (ускорение свободного падения g). Введем угол отклонения маятника от вертикали. Движение происходит по дуге окружности, скорость тела v = l. Уравнение динамики вращательного движения относительно точки подвеса ml2 2 = mgl sin. Здесь ml2 момент инерции маятника. Рассматривая угол и скорость как функции времени, получаем систему из двух уравнений:

При малых углах отклонения можно считать sin, и мы получаем уравнение гармонических колебаний d2 /dt2 + 0 = 0. Частота 0 = g/l в этом случае не зависит от размаха колебаний. В общем случае, когда угол отклонения не мал, эту систему дифференциальных уравнений приходится решать на компьютере.

1.3.1. Безразмерные переменные В численном эксперименте следует избегать появления очень больших и очень малых чисел. Поэтому вовсе не обязательно привязываться к какой-либо стандартной системе единиц. Наоборот, удобно считать основные параметры задачи единичными величинами. Таких в динамических задачах можно выбрать три (так как имеется три основные величины: масса, длина и время). В данном случае удобно положить m = 1, l = 1, g = 1. Для того, чтобы перевести значения остальных переменных (кроме основных) в размерные величины, необходимо выразить их размерности через размерности основных и затем подставить численные значения основных величин. Например, скорость будет измеряться в единицах v0 = gl, энергия в единицах E0 = mgl, время в единицах t0 = l/g и т. д.

1.3.2. Фазовые диаграммы В процессе движения маятника угол отклонения и скорость меняются со временем. Дополнительную информацию можно получить, изображая состояние системы в координатах угол – момент импульса (, L).

Здесь L = ml2 d/dt. Такая система координат называется фазовой плоскостью, а линия движения состояния системы фазовой траекторией. Углы отклонения маятника и + 2 физически неразличимы, поэтому можно рассматривать угол в конечных пределах (например, < ). При этом надо представлять себе, что точки фазовой плоскости (, L) и (, L) отождествлены, то есть края фазовой плоскости = ± склеены между собой.

При малых амплитудах колебания являются гармоническими /max = cos(t + ), v/vmax = sin(t + ), фазовая траектория представляет собой окружность в координатах (/max, v/vmax ). Точка, показывающая текущее состояние маятника, движется по этой окружности с постоянной скоростью. Вообще, на фазовой плоскости любым периодическим движениям соответствуют замкнутые кривые.

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

1.3.3. * Физический маятник Физический маятник представляет собой тело конечных размеров, подвешенное в некоторой точке. Разобьем мысленно тело на маленькие кусочки, имеющие массы mi и находящиеся в точках ri, считая от точки подвеса.

Здесь и далее полужирным шрифтом будем обозначать векторы. Тогда полная масса тела m = mi, радиус-вектор центра масс R = mi ri /m, момент инерции тела относительно горизонтальной оси вращения, проходящей через точку подвеса, I = mi h2, где hi перпендикуляр от точки mi до оси вращения. В положении равновесия вектор R направлен вертикально вниз (вдоль g). Уравнение движения для физического маятника где M = mgR sin момент силы тяжести относительно точки подвеса. Например, для однородного стержня массой m и длиной l, подвешенного за один из концов, I = ml2 /3, R = l/2 и уравнение движения d2 /dt2 = 3g sin /2l. В общем случае это уравнение имеет тот же вид, что и для математического маятника, только частота колебаний другая.

Например, при малых колебаниях стержня 0 = 3g/2l.

Задания 1. Выясните, как меняется вид фазовой траектории при увеличении энергии колебаний. Найдите сепаратрису. Постройте график зависимости периода колебаний от амплитуды.

2. На реальный маятник действует сила трения. При не слишком больших скоростях можно принять Fтр = kv. Что происходит при этом с фазовой траекторией? При каких значениях коэффициента k колебания не возникают?

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

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

1.4. Движение в центральном поле (планета) Программа planet.pas.

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

Рассмотрим планету массой m, движущуюся вокруг звезды массой M. Сила притяжения по модулю равна где r = rm rM вектор расстояния между телами, и направлена вдоль r. Считаем M m, тогда ускорение звезды очень мало, и ее движением можно пренебречь. Направление начальной скорости планеты и линия, соединяющая планету со звездой, образуют плоскость, в которой происходит все дальнейшее движение. Назовем ее плоскостью XY и поместим звезду в начало координат. Радиус-вектор планеты r = (x, y).

Вместо произведения GM можно ввести параметр = GM. Переходя к безразмерным переменным, удобно считать = 1. Сила, действующая сы планеты. Поскольку сила всегда направлена по радиусу, ее момент M = [r F] равен нулю, и момент импульса L = [r p] сохраняется (второй закон Кеплера). Движение планет происходит по эллипсам (первый закон Кеплера), а квадраты периодов обращения относятся как кубы больших полуосей (третий закон Кеплера). Потенциальная энергия планеты равна (полагая ее равной нулю на бесконечном расстоянии) Первая космическая скорость vI соответствует вращению по круговой орбите. Вторая космическая скорость vII минимальная скорость, необходимая, чтобы улететь на бесконечность. Хорошо известно, что для тела в поле тяготения Земли вблизи поверхности vI = 7,9 км/c, vII = 11,2 км/с. (Иногда говорят еще про третью космическую скорость, это понятие возникает при рассмотрении движении тела одновременно под влиянием Земли и Солнца. Третья космическая скорость минимальная скорость, необходимая для того, чтобы, взлетев с Земли, навсегда покинуть Солнечную систему vIII = 16,4 км/с. ) При рассмотрении движения планеты в поле тяготения Солнца можно обобщить понятие первой космической скорости vI и получить ее значение, приравняв ускорение планеты за счет силы притяжения Солнца центростремительному ускорению при ее движении по круговой орбите радиуса r:

Отсюда vI = GM/r. Обратите внимание на то, что первая космическая скорость зависит от положения тела относительно центра тяготения.

Чтобы найти vII, заметим, что полная энергия планеты E = mv 2 / GM m/r сохраняется. На бесконечности должно быть E = mv 0, равенство как раз соответствует второй космической скорости для планеты (эту скорость еще называют скоростью отрыва). Таким образом, и vII = 2GM/r. Видно, что vII = 2vI.

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

Сила, действующая между двумя электрическими зарядами q1, q2, находящимися на расстоянии r, определяется законом Кулона F = q1 q2 /r2 и имеет тот же вид, что и гравитационное взаимодействие. Рассматривая в качестве тяжелой частицы протон, а в качестве легкой электрон, получим так называемый ридберговский атом. В таком атоме расстояние между протоном и электроном достаточно большое, и электрон можно считать классической частицей. Очень интересно изучить поведение таких атомов в электрическом и магнитном полях.

Задания 1. Определите методом подбора первую космическую скорость. Запишите ее значение. Полезно сделать так, чтобы при нажатии определенной клавиши планета устанавливалась в начальное положение, имея нулевую начальную скорость по оси X, а вы попадали бы сразу в таблицу для того, чтобы изменить значение начальной скорости по оси Y. Разумно проверять значение переменной key, которой присваивается код последней нажатой клавиши после выполнения оператора key:=CheckPoint(), затем использовать процедуры DrawTable(MainTable) для отрисовки таблицы и ModifyTable(MainTable) для ее изменения.

2. Подберите вторую космическую скорость. Изобразите графики кинетической, потенциальной и полной энергии в зависимости от r.

Во сколько раз vII больше vI ? Что происходит с графиком полной энергии при увеличении шага по времени?

3. Вычислите значение момента импульса, проверьте его сохранение, выведите график. Абсолютное значение L = rp sin, угол между векторами r и p. Как выражается L через компоненты 4. Проверьте выполнение третьего закона Кеплера.

5. Рассмотрите движение частицы в поле вида F rn. При каких n получаются замкнутые орбиты? В чем их сходство и отличие?

6. Учтите влияние на спутник торможения в атмосфере (простейший способ ввести силу трения Fтр = kvv). Почему торможение приводит к увеличению кинетической энергии?

7. Рассмотрите движение одной частицы при наложении постоянной силы вдоль оси Z (это может быть, к примеру, электрическое поле).

Дополнительная сила должна быть мала по сравнению с гравитационной, Fэл Fгр. Выведите на экран вид орбиты в трех проекциях XY, XZ, Y Z в разных окнах. Объясните, что происходит с орбитой. К чему приводит увеличение силы?

8. Как движется частица, если направление постоянной силы лежит в плоскости орбиты (классический эффект Штарка)? Желающим более подробно разобраться рекомендуем посмотреть [7], задачи 2.35, 9. Рассмотрите движение частицы при добавлении магнитного поля, направленного вдоль оси X (или Z) (классический эффект Зеемана). При движении заряженной частицы в магнитном поле возникает сила Лоренца Fл = q[v B]/c. Как выражаются компоненты 10.* Существует теорема Лармора о том, что влияние слабого магнитного поля (магнитная сила гораздо меньше силы притяжения) эквивалентно переходу в систему отсчета, вращающуюся с Ларморовской частотой =. Проверьте это утверждение.

11. Рассмотрите обратный случай магнитная сила гораздо больше силы притяжения (или Ларморовская частота много больше невозмущенной частоты обращения, 2/T ). Как выглядит траектория частицы в этом случае?

12. Добавьте еще одну планету, не взаимодействующую с первой. Как выглядит движение второй планеты при наблюдении ее с первой планеты? Чем отличаются картины, наблюдаемые с внутренней и внешней планет?

13. Введите взаимодействие между планетами. Подберите параметры так, чтобы получилась устойчивая система типа Земля – Луна, которая движется вокруг Солнца.

* Примечание. Для двух векторов a = (a, a, a ) и b = (b, b, b ) их векторное произведение c = [a b] = (cx, cy, cz ) имеет компоненты 1.5. Случайные блуждания Программа dius.pas.

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

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

Хаотическое движение маленьких частиц вещества можно наблюдать в микроскоп. Первым такое наблюдение сделал Р. Броун (R. Brown) в 1827 году. В его честь такое движение называется броуновским (правильнее было бы брауновским). Полная теория броуновского движения дана в 1905–1906 гг. А. Эйнштейном и М. Смолуховским.

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

1.5.1. Одномерные случайные блуждания Рассмотрим простую модель пусть частицы могут двигаться только вдоль прямой, делая шаги случайной длины. Примем также, что величина и направление каждого шага определяются независимо от положения частицы и от предыдущих шагов (модель пьяного моряка ). Будем наблюдать за частицей через равные промежутки времени. Координата частицы вычисляется по рекуррентной формуле xk = xk1 + k, где k очередной шаг блуждания. При отсутствии силового поля смещение влево и вправо равновероятны. В общем случае вероятность того, что длина шага лежит в промежутке от до + d, равна dp = w()d. Функция w называется плотностью вероятности для величины шага, который называют случайной величиной. Примем вначале простейшую зависимость симметричное равномерное распределение:

В языке Pascal имеется стандартный генератор случайных чисел random, который дает случайное вещественное число, равномерно распределенное в интервале от 0 до 1. С его помощью такую плотность вероятности для можно получить, если вычислять = 2 · random 1.

Очевидно, что полная вероятность сделать хоть какой-нибудь шаг равна Среднее смещение за шаг по времени = w()d = 0 следует из того, что шаги влево и вправо равноправны. Нам понадобится еще средний квадрат смещения за один шаг. Так как 1 < < 1, то 2 = 2 w()d = 1/3. (Вообще, среднее значение любой функции f () вычисляется как f () = f ()w()d.) 1.5.2. Средние значения Рассмотрим положение частицы в какой-либо момент времени. Пусть частица начинает движение в точке x0 = 0 и после k шагов имеет координату xk = xk1 + k. Ее квадрат x2 = x2 + 2xk1 k + k. Мы можем повторить опыт с частицей много раз (или выпустить сразу много частиц) и найти среднее значение координаты xk, которое равно xk = = xk1 + k, так как среднее значение суммы нескольких величин равно сумме их средних значений. Поскольку k = = 0, то xk = = x0 = x0 = 0.

Теперь найдем средний квадрат координаты x2 = x2 + + 2 xk1 k + k. Случайные величины xk1 и k независимы, тогда среднее значение их произведения равно произведению средних xk1 k. Поэтому второе слагаемое исчезает, остается x2 = x2 + + k = k 2, так как x2 = 0. Номер шага пропорционален времени k t, отсюда получаем известную формулу диффузии x2 = 2Dt (в одномерном случае). В нашей простейшей модели блужданий коэффициент диффузии получается равным D = 2 /2 = 1/6.

1.5.3. Функция распределения Рассмотрим ансамбль частиц, которые одновременно начинают совершать одномерные случайные блуждания из начального положения x = 0.

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

Если мы будем одновременно увеличивать число частиц и число столбцов гистограммы (при этом в каждый столбец гистограммы должно попадать много частиц), то в пределе гистограмма приближается к гладкой кривой, называемой функцией распределения f (x). Вероятность найти частицу на малом интервале от x до x + dx пропорциональна произведению f (x)dx. Например, после первого шага точная функция распределения наших частиц совпадает с w(x).

После k шагов координата частицы xk = i. В теории вероятноi= стей есть так называемая центральная предельная теорема: распределение случайной величины, значение которой является суммой большого числа случайных величин, стремится к распределению, называемому нормальным. Более подробно об этом можно прочитать в [4], гл. 11.

В нашем случае, для процесса диффузии частиц функция распределения на k-шаге по времени принимает вид нормального распределения:

если все частицы в момент времени t = 0 находились в начале координат. Строго говоря, эта функция близка к полученному распределению 1.5.4. Многомерные случайные блуждания Предоставим теперь нашим частицам возможность двигаться также по координате y рассмотрим двумерные случайные блуждания (случай трех измерений получается аналогично). Можно независимо задавать смещение по вертикали y равномерно распределенным таким же образом, как и смещение по горизонтали x.

Функция распределения в двумерном случае представляется в виде произведения двух функций распределения по координатам x и y, так как x и y являются независимыми случайными величинами. Удобно перейти к распределению по радиусу r = x2 + y 2. Так как произведение плотностей вероятности f (x) и f (y) зависит только от r то вероятность найти частицу, расстояние которой от начала координат лежит в интервале от r до r+dr (тонкое кольцо площадью dS = 2rdr), равна Видно, что при малых r эта вероятность возрастает линейно от нуля, потому что экспоненциальный множитель мало отличается от единицы.

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

В трехмерном случае вместо кольца надо рассмотреть тонкую сферу, объем которой dV = 4r2 dr. В результате получаем При использовании равномерных распределений для x и y получается некоторая анизотропия, особенно сначала (после первого шага частицы заполняют квадрат, а не круг), однако она быстро сглаживается. Избавиться от анизотропии можно, если задавать x и y независимо равномерно распределенными от до, а затем отбрасывать точки, оказавшиеся вне круга с радиусом. Другой способ задание равномерно распределенных случайных величин длины шага r и угла = 2 · random приводит к неравномерному распределению, плотность числа частиц уменьшается при увеличении r. Равномерную плотность можно получить, если задавать шаг по радиусу в виде r = random (подробности в [4], глава 10).

1.5.5. Случайные блуждания во внешнем поле.

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

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

Такая модель соответствует распределению плотности в изотермической атмосфере. Частица массой m имеет в поле тяжести на высоте h потенциальную энергию U = mgh. В равновесии концентрация подчиняется распределению Больцмана где k = 1, 38·1016 эрг/К постоянная Больцмана. Характерная энергия теплового движения kT и средний квадрат случайного смещения пропорциональны. Найдем среднюю высоту облака частиц Для атмосферы Земли h составляет примерно 8 км эта величина мала по сравнению с радиусом Земли, поэтому можно считать g постоянным. Реально плотность падает быстрее из-за уменьшения температуры с высотой. Интересно заметить, что в изотермической атмосфере произведение mg h = kT не зависит от g.

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

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

3. Постройте одновременно с гистограммой график теоретической функции распределения (процедура GraFun). Согласуйте масштаб для функции с масштабом гистограммы.

4. Реализуйте в программе двумерные случайные блуждания, в нулевой момент времени все частицы расположены в точке (0,0). Постройте график среднего квадрата радиуса, среднего квадрата x от времени. Постройте гистограмму распределения частиц по радиусу.

Как объясняется спад к нулю при малых r (в области, где частицы расположены наиболее плотно)?

5. Рассмотрите одномерные случайные блуждания в поле тяжести с упругим отражением частиц от нижней стенки. Постройте гистограмму, определите среднюю высоту облака частиц h в зависимости от g. Как ведет себя выражение h g? Для ускорения работы можно ввести два сорта частиц с разными массами (разные g) и вывести одновременно их распределения (для вывода вторых графиков используйте процедуры TraceB и HistDrawT).

6.* Другой, более подробный, способ описания броуновского движения уравнение Ланжевена:

На частицу действует вязкое трение и случайная сила, возникающая в результате ударов молекул окружающей среды. За время t действие сил приводит к изменению скорости Среднее значение случайной силы F = 0, ее средний квадрат связан с коэффициентами вязкости и диффузии F2 = 3dD 2 /t.

Здесь d размерность пространства, t шаг по времени. (Более правильно говорить о корреляционной функции силы, подробности можно посмотреть в [3].) Исследуйте, как меняется со временем средняя координата частицы, средний квадрат координаты.

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

Глава 2.

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

1. Неравновесная агрегация, фрактальные кластеры.

2. Электрический пробой.

3. Молекулярная динамика.

4. Решеточные газы, решеточное уравнение Больцмана.

5. Химические реакции, горение (стохастическое).

6. Теплопроводность, детерминированное горение.

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

8. Колебания цепочек.

9. Рост дендритов.

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

2.1.1. Агрегация, ограниченная диффузией Простейший случай агрегация, ограниченная диффузией (Diusion Limited Aggregation, DLA [8]). Для начала рассмотрим сеточную модель.

Возьмем регулярную сетку на плоскости, например, квадратную.

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

Некоторые указания. Для ускорения работы программы разумно выпускать частицы с круга радиусом немного больше Rmax текущего максимального радиуса агрегата. Угловую координату естественно задавать равномерно распределенной на интервале от 0 до 2, используя формулу Если частица уходит далеко (например, дальше, чем 2Rmax ), уничтожаем ее и выпускаем новую.

2.1.2. Фрактальная размерность Фигура на плоскости или тело в пространстве имеют размерность. Определить ее можно разными способами. В случае, когда у фигуры есть выделенная центральная точка, можно построить много сфер различного радиуса с центром в этой точке. Для каждой сферы можно вычислить массу части фигуры, которая оказалась внутри этой сферы. В случае, когда масса пропорциональна радиусу сферы в некоторой степени (m RD ), показатель степени D называется размерностью тела. Это так называемый метод сфер или ящиков. Для линий D = 1, у плоских фигур D = 2, у обычных тел D = 3. Однако, многие объекты в природе имеют размерность, выражающуюся дробным числом (например, у кластера DLA на плоскости размерность D 1,71 ± 0,02, в трехмерном пространстве D 2,50). Такие тела, следуя Б. Мандельброту [9], называют фракталами (от латинского слова fractus дробный). Фракталами являются также дендриты, вырастающие при электроосаждении металлов, кластеры, полученные при агрегации коллоидов; фрактальную структуру имеют ветви деревьев, кровеносная система и т. п.

Удобно строить график ln m в зависимости от ln R. Если он близок к прямой, тангенс ее наклона как раз определяет D. Средняя плотность таких агрегатов a m/Rd RDd убывает с увеличением радиуса (здесь d размерность пространства, D < d).

Другой способ определения размерности метод подсчета клеток (box counting). Возьмем куб (шар) со стороной L0, который целиком заключает в себе объект. Разобьем куб на меньшие, со стороной L1 (к примеру, хотя и не обязательно, L1 = L0 /2). Пусть число непустых таких кубов (в которые попадает часть объекта) будет N1. Продолжим процедуру для L2 и N2. Если оказывается, что Ni LD, то размерность объекта равна D. К примеру, размерность береговой линии Норвегии 1,52 ± 0,01 [10]. Размерности, определенные этим и предыдущим способом, совпадают для простых тел и могут слегка различаться для хитрых кластеров.

Еще один способ может применяться при наблюдении за процессом роста агрегата от центра. В этом случае число частиц в кластере N Rch. В качестве характерного радиуса Rch можно выбирать, например, максимальный радиус кластера Rmax, однако более правильно использовать радиус гирации Rg = рассмотрении момента инерции кластера, состоящего из частиц единичной массы, N Rg = ri.

* Экспериментально можно определять размерность агрегатов методом малоуглового рассеяния света, рентгеновских лучей или нейтронов. Зависимость интенсивности рассеяния от угла определяется выражением где s длина вектора рассеяния, угол рассеяния, длина волны излучения. Эта формула справедлива при s a, где a размер частиц, из которых состоит кластер. Для бльших s начинает преобладать рассеяние отдельными частицами, и I(s) s4.

Рис. 2.1. Множество Кантора (а), кривая Коха (б), треугольник Серпинского (в) Примеры математических фракталов 1. Множество Кантора. Возьмем отрезок и выкинем среднюю его треть. Далее у двух оставшихся частей также удаляем среднюю треть. Так продолжаем до бесконечности. На рис. 2.1,а показаны первые пять итераций. То, что остается, называется множеством Кантора ( Канторовой пылью ). Его размерность легко определяется методом подсчета клеток. Удобно взять Li+1 = Li /3, тогда 2. Кривая Коха. Снова начинаем с отрезка, среднюю треть которого заменяем на два кусочка такой же длины, но расположенных под углом 60. Поступаем так же с каждым из получившихся отрезков. Размерность получившейся снежинки ln 4/ ln 3 1, (рис. 2.1,б).

3. Треугольник Серпинского. Берем равносторонний треугольник и соединяем середины сторон. Внутренний маленький треугольник выкидываем. Продолжаем процесс. Размерность оставшейся салфетки ln 3/ ln 2 1,585 (рис. 2.1,в).

4. Траектория броуновской частицы. При случайных блужданиях смещение частицы за время t изменяется по закону r t (см.

раздел 1.5). При этом длина пути L = va t, где va средняя скорость движения. То есть, L r, D = 2, если не учитывать самопересечения.

Желающим более подробно познакомиться с фракталами рекомендуем книгу [10].

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

2.1.4. Химически-ограниченная агрегация При диффузионно-ограниченной агрегации частица всегда прилипает к кластеру с вероятностью 1. Можно уменьшить вероятность прилипания.

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

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

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

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

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

Задания 1. Напишите программу модели DLA на сетке.

2. Постройте график зависимости ln N (числа частиц в кластере) от ln Rg (радиуса гирации). Определите размерность агрегата.

3. Определите размерность кластера методом ящиков.

4. Найдите зависимость размерности от вероятности прилипания частицы. Указания. Генератор случайных чисел random выдает числа, равномерно распределенные между 0 и 1. Поэтому условие random < p будет выполняться с вероятностью p, если p 1.

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

5. Рассмотрите случай химически-ограниченной агрегации. Можно задать, например, следующие вероятности прилипания в зависимости от числа соседей: p1 = 0,01, p2 = 0,1, p3 = 0,9. Определите размерность полученного кластера.

6. Реализуйте бессеточную модель DLA.

7. Напишите программу моделирования баллистической агрегации.

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

8. Рассмотрите кластер-кластерную агрегацию. Какова размерность получившегося агрегата?

2.2. Электрический пробой Человек познакомился с искровым разрядом задолго до того, как приступил к научному познанию мира. Феерическое и грозное явление природы молния с точки зрения физики являет собой пример грандиозного искрового разряда в атмосфере. Началом систематического исследования электрического разряда можно считать экспериментальные исследования Б. Франклина в середине XVIII века, в которых было доказано единство природы молнии и лабораторной электрической искры.

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

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

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

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

При напряжениях, более высоких, чем те, которые приводят к образованию короны, в газах возникают так называемые стримеры. Стримеры представляют собой систему слабосветящихся проводящих каналов, образующуюся в газе в области наиболее сильного электрического поля. Стример прорастает, как правило, с одного из электродов и при высоких напряжениях может ветвиться. Ветвление стримера происходит нерегулярно, и на сегодняшний день можно считать, что эти ветвления носят случайных характер. Разветвленную стримерную вспышку часто называют импульсной короной. Размер области, в которой развиваются стримеры, может составлять несколько метров даже в лабораторных условиях. Скорость продвижения стримера вглубь межэлектродного промежутка не меньше 10 км/с и может достигать 10000 км/с. Вспышкой импульсной короны начинается искровой разряд в воздухе или других газах в длинных промежутках между электродами. В случае пробоя в газах при определенных условиях энерговыделения в стримерных каналах они превращаются в так называемые лидерные каналы. Лидер это плазменное образование очень высокой светимости и настолько высокой проводимости, что его в некотором смысле можно считать продолжением электрода. Скорость распространения лидера по порядку величины составляет 10 км/с. Перед головной частью лидера образуется стримерная корона, от которой зависит дальнейшая динамика лидерного канала.

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

На сегодняшний день менее всего изучен механизм роста и ветвления стримеров при электрическом разряде. Согласно современным представлениям, рост кончика стримера определяется величиной напряженности электрического поля перед ним. Величина электрического поля зависит не только от падения напряжения между электродами, но и от радиуса стримера и от скорости его роста. Эти три величины радиус, локальная напряженность электрического поля и скорость роста кончика стримеГлава 2. Задания для самостоятельной работы ра связаны между собой. При определенных значениях этих величин кончик стримера может разветвляться. В 1984 г. была предложена первая модель стохастического роста стримеров, в которой рост их ветвей связывался с локальным электрическим полем [11]. Основные положения модели изложены ниже. Модель легко реализуется на компьютере и позволяет качественно воссоздать картину системы стримеров, наблюдаемых при искровом разряде в газах. Авторы модели Нимейер, Пьетронеро и Висман (НПВ) использовали ее для описания пространственной структуры газового разряда, возникающего на поверхности диэлектрика, помещенного между двумя проводниками, в газе. Картины такого поверхностного разряда наблюдались впервые в 1777 г. и получили название фигур Лихтенберга. Предлагается смоделировать возникновение фигур Лихтенберга с помощью моделей стохастического роста стримеров.

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

По теореме Гаусса поток вектора индукции электрического поля D через любую замкнутую поверхность S равен нулю при отсутствии внутри поверхности свободных электрических зарядов (n – вектор внешней нормали к поверхности). Для большинства диэлектриков индукция электрического поля выражается через электрическое поле как D = E. В общем случае вектор E имеет три компоненты (Ex, Ey, Ez ).

Рассмотрим в пространстве кубическую решетку с ячейками со сторонами h по всем координатам x = h, y = h и z = h. Сначала рассмотрим только один ряд ячеек вдоль оси x (рис. 2.2). Пусть электрический потенциал принимает в центре i-той ячейки значение i,j,k.

Тогда проекция электрического поля на ось x на левой грани ячейки приближенно равна Ex (x) (i,j,k i1,j,k )/h, а на правой грани Ex (x + x) (i+1,j,k i,j,k )/h. Заметим, что эти формулы являются конечно-разностными аппроксимациями первых производных. Поток вектора индукции электрического поля изнутри ячейки через эти две грани равен x = Ex (x + x)h2 Ex (x)h2 = (i+1,j,k 2i,j,k + i1,j,k )h.

Аналогично поток через нижнюю и верхнюю грани равен y = Ey (y + y)h2 Ey (y)h2 = (i,j+1,k 2i,j,k + i,j1,k )h, а через переднюю и заднюю грани z = Ez (z + z)h2 Ez (z)h2 = (i,j,k+1 2i,j,k + i,j,k1 )h.

Складывая эти уравнения, вычислим полный поток изнутри ячейки и поделим его на объем ячейки (i,j+1,k 2i,j,k + i,j1,k ) (i,j,k+1 2i,j,k + i,j,k1 ) Каждая дробь представляет собой разностную аппроксимацию второй частной производной от потенциала по соответствующей координате.

Так как в данном случае полный поток равен нулю, то при h 0 уравнение для потенциала в однородном диэлектрике сводится к хорошо известному уравнению Лапласа:

Для краткости уравнение Лапласа часто записывают в символической форме 2 = 0, где оператор 2 = 2 + 2 + 2.

Используя условие равенства нулю полного потока из уравнения (2.1) можно также получить уравнение i1,j,k + i+1,j,k + i,j1,k + i,j+1,k + i,j,k1 + i,j,k+1 6i,j = 0. (2.2) Далее будем рассматривать плоский случай. Потенциал изменяется только в плоскости XY, поэтому по теореме Гаусса для квадратной ячейки сетки с номером i, j вместо формулы (2.2) получим Последнее уравнение можно переписать в виде:

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

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

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

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

Тогда потенциал второго электрода равен приложенному напряжению.

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

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

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

Зная потенциал, можно вычислить электрическое поле. Среднее значение проекции электрического поля на звено, соединяющее узлы A и B, EAB = (A B )/lAB. Нам необходимо знать поле на звеньях, выходящих из стримерной структуры. Считая ее потенциал нулевым, получаем просто |E| = B для горизонтальных и вертикальных звеньев, |E| = B / 2 для диагональных. Здесь индекс B означает узлы, соседние со структурой.

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

Модель НПВ. Нимейером, Пьетронеро и Висманом [11] впервые была предложена модель, которая позволяет описать рост структур разряда в диэлектриках. В основе модели лежит предположение, что структура растет случайным образом, причем вероятность роста зависит только от локального электрического поля вблизи структуры. Эта модель достаточно легко реализуется на компьютере. Обычно моделируется разряд в диэлектрике, помещенном между двумя электродами, разность потенциалов между которыми V.

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

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

На каждом шаге роста случайный процесс выбора новой веточки структуры реализуется следующим алгоритмом. Пробегаем по всем M узлам решетки, в которые возможен рост, и рассчитываем сумму где величина Ek своя для каждой пары узлов и находится на каждом шаге роста по текущему распределению электрического поля. Затем разыгрывается случайное число, равномерно распределенное от 0 до Z. Для этого используется формула = Z ·random. Затем повторно шаг за шагом рассчитывается (2.4) до тех пор, пока текущая сумма S Ek не станет больше. Тот узел, для которого сумма стала больше, присоединяется к структуре. Новой образовавшейся веточке присваивается значение потенциала того электрода, с которого начался рост этого стримера. Таким образом формируется эквипотенциальная структура.

В этой модели каждый шаг роста во времени имеет номер, но, к сожалению, номеру не сопоставлен реальный интервал времени t. Структура разряда растет до тех пор, пока не достигнет противоположного электрода.

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

Флуктуационный критерий роста. В многозвенных моделях, наоборот, пренебрегается влиянием друг на друга проводящих звеньев, возникающих на данном шаге по времени. Пусть пробой происходит в тех областях диэлектрика, где величина поля превосходит некоторое пороговое значение E (электрическую прочность). Однако это не совсем правильно, так как реальный диэлектрик не является полностью однородным. Существуют слабые места, связанные, например, с повышенной ионизацией за счет космических частиц, с флуктуациями температуры и т. д. То есть, критерий возникновения пробоя можно представить в виде E > E, где величина флуктуации [12, 13]. Очевидно, что большие значения будут встречаться реже, чем маленькие. Для такого случая независимых редких событий распределение для флуктуаций экспоненциальное f () exp(/g)/g, где g некоторый параметр.

Такое распределение вероятностей можно получить, задавая случайное число = g ln(random) (подробности в [4], глава 10). Это соответствует вероятности пробоя p(E) = exp((E E )/g).

Развитие пробоя происходит следующим образом: проверяются все звенья, выходящие из стримерной структуры. Те из них, для которых Ei > E i, присоединяются к структуре. Эта модель роста принадлежит к классу многозвенных. Шаг роста во времени t считается постоянным.

Модели стохастического времени запаздывания. В [14] для каждого возможного звена было введено случайное время ожидания пробоя i = ln(random)/r(E). Случайные числа отражают стохастический (вероятностный) характер процесса. Тогда вероятность пробоя за малый интервал времени 0 равна p(E) = 0 r(E). Функция вероятности роста r(E) может быть любой быстро растущей функцией от локального электрического поля. В однозвенной модели пробивается звено, у которого время i минимально [14]. Шаг роста во времени принимается равным этому минимальному значению i.

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

Как уже говорилось, вероятность пробоя резко возрастает с увеличением напряженности электрического поля. Поэтому часто используется степенная аппроксимация p(E) E. При = 2 вероятность пробоя пропорциональна мощности локального энерговыделения w = E 2, где характерная величина электропроводности плазмы в каналах структуры.

Интересно отметить, что такой вид p(E) может описывать не только электрический пробой, но и другие процессы. Например, при = вероятность роста не зависит от поля, таким образом можно описать рост колоний бактерий. При = 1 вероятность p, что соответствует диффузионно-ограниченному росту (модель DLA, см. раздел 2.1).

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

Задания 1. Напишите программу вычисления электрического потенциала итерационным методом.

2. Рассмотрите пробой в геометрии острие – плоскость с использованием флуктуационного критерия роста.

3. То же в геометрии точка – окружность. Как меняется густота ветвей с радиусом стримерной структуры? Для этого можно рассмотреть отношение числа точек пересечения структуры с окружностью некоторого радиуса к длине этой окружности.

4. Реализуйте в геометрии острие – плоскость или точка – окружность однозвенную или многозвенную модель со степенной зависимостью вероятности роста от напряженности поля p E. Рассмотрите случаи = 0, 1, 2. Как меняется геометрия стримерной 2.3. Молекулярная динамика Метод молекулярной динамики (МД) рассматривает поведение вещества на микроуровне мы наблюдаем за движением отдельных молекул. При этом мы хотим понять поведение сложной многочастичной системы. Ясно, что макроскопические тела содержат слишком много молекул, чтобы их можно было смоделировать. Однако, применение метода МД даже к небольшим системам, состоящим из нескольких сотен или тысяч частиц, дает много для понимания наблюдаемых свойств газов, жидкостей и твердых тел. Кроме того, малые системы (наночастицы) имеют свои совершенно новые свойства, изучение которых очень интересно.

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

2.3.1. Уравнения МД Несмотря на то, что микрочастицы (атомы и молекулы) квантовые объекты, в некоторых практически интересных случаях возможно рассматривать их классически. Тогда влияние молекул друг на друга описывается потенциалом взаимодействия U. Простейший случай парное взаимодействие, когда Uij = U (rij ), где rij = |ri rj |. Суммарная потенциальная энергия равна Ep = Uij. Сила, действующая на данную чтобы он обращался в ноль на бесконечности U (r) = ULJ (r) ULJ (rc ).

2.3.3. Вычисление ускорений Как уже было сказано, расчет взаимодействий занимает основное время в ходе вычислений. Поэтому оптимизировать нужно именно его. Вдвое увеличить скорость можно, учитывая третий закон Ньютона Fij = Fji.

Примерный вид куска программы на языке Pascal имеет вид:

for i := 1 to N-1 do for j := i+1 to N do begin {Вычисление Fij} a[i] := a[i]+Fij/m[i];

a[j] := a[j]-Fij/m[j];

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

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

Несмотря на непрерывно возрастающую мощность компьютеров, все же число частиц, доступное моделированию, гораздо меньше характерного числа молекул в макроскопическом объеме ( 1022 ). Чтобы соотносить результаты расчета со свойствами гораздо большей системы, необходимо исключить влияние границ расчетной ячейки, которое пропорционально N 1/2 в плоском случае и N 1/3 для трех измерений. Простым способом является задание периодических граничных условий. При этом края расчетной области замыкаются в тор, т. е. частица, вылетающая за левую границу, появляется справа, вылетающая вверх снизу и т. п.

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

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

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

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

Например, для плоского кластера получаем:

где vcm скорость центра масс кластера.

Давление можно вычислять различными способами. Один из них следует из теоремы вириала:

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

Другая величина, которая может служить индикатором начала плавления кластера среднеквадратичная флуктуация длины связи. Она вычисляется следующим образом течение по крайней мере нескольких шагов частица будет взаимодействовать только с этими соседями. Этот метод называется методом списков Верле. При составлении списков разумно пользоваться ячейками. В литературе рекомендуется брать rv = 1,05rc и обновлять списки всякий раз, когда смещение какой-либо частицы превысит 0,5(rv rc ) [21]. В трехмерном случае использование списков Верле сокращает число расчетов раз в 6, в двумерном варианте ускорение не столь значительно.

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

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

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

Задания 1. Напишите программу двумерной молекулярной динамики. Задайте начальные условия с нулевым суммарным импульсом частиц. Граничные условия задайте периодическими. Подберите подходящий шаг по времени. Проверьте сохранение полной энергии и импульса.

2. Постройте гистограмму распределения частиц по скоростям. Сравните ее с распределением Максвелла.

3. Исследуйте в двумерной модели плавление и затвердевание малых кластеров с магическим числом атомов в них 7, 19, 37, 61 и т. д.

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

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

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

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

5.* При нескольких значениях плотности найдите зависимости P от T. По полученным зависимостям оцените постоянные a и b в уравнении состояния Ван-дер-Ваальса. Определите температуру, плотность и давление в критической точке.

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

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

2.4. Решеточные газы, решеточное уравнение Больцмана 2.4.1. Решеточные газы (Lattice-Gas Automata, LGA) Рассмотрим квадратную решетку, в узлах которой могут обитать частицы единичной массы. Выберем расстояние между узлами x в качестве единицы длины, а шаг по времени t за единицу времени. Каждая из частиц может иметь скорость, направленную в один из соседних узлов, такой величины, что за один шаг по времени частицы могут перелететь как раз в соседний узел (рис. 2.6,а). В узле может быть не более одной частицы с данным направлением скорости (принцип исключения).

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

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

Благодаря принципу исключения наличие частицы, имеющей скорость по каждому направлению, может быть закодировано одним битом (0 нет частицы, 1 есть). Так можно записать состояние каждого узла в четырех битах. Можно пронумеровать направления в порядке возрастания угла. Например, в двоичной системе первое состояние на рис. 2.6,б записывается как S = 01012. Всего различных состояний в узле может быть 16, включая S = 0 (нет ни одной частицы). Все операции сводятся к целочисленной арифметике, это означает высокую скорость расчетов и отсутствие ошибок округления. Кроме того, все вычисления локальные, поэтому их можно выполнять параллельно.

Обозначим возможные направления скорости d1, d2, d3, d4 (рис. 2.6,а).

Рис. 2.6. Возможные направления скорости частиц в модели HPP (а) и возможные столкновения, в которых скорости частиц изменяются (б) Рис. 2.7. Решетка и некоторые возможные столкновения частиц в модели FHP-I (а), некоторые возможные столкновения с участием покоящихся частиц в модели FHPIII (б) Тогда d1 = 00012 = 1, d2 = 00102 = 2, d3 = 01002 = 4, d4 = 10002 = 8.

Приведем в явном виде выражения для основных операций.

1) Добавление к состоянию S частицы с направлением скорости dk :

2) Проверка: есть ли в состоянии S частица с направлением скорости if (S and dk ) = 0 такая частица есть, иначе нет.

Здесь or двоичная побитовая операция или, а and двоичная операция и.

2.4. Решеточные газы, решеточное уравнение Больцмана Недостатки модели HPP связаны с тем, что квадратная сетка с 4 возможными направлениями скорости частиц недостаточно симметрична.

Однако, уже простейшая модель LGA на треугольной сетке с 6 возможными векторами скорости частиц в узле (модель FHP-I, [24]) обладает достаточной симметрией. В более сложных вариантах возможно наличие покоящихся частиц (модель FHP-III, [24]). Геометрия решетки и возможные столкновения частиц для этих моделей показаны на рис. 2.7.

Другой способ устранить недостатки квадратной сетки ввести еще возможность движения частиц по диагоналям (скорость 2). Вместе с покоящимися частицами получаем 9 направлений скорости [25].

Так как модули скоростей различны, возможен нетривиальный закон Рис. 2.8. Геометрия и примеры столкносохранения энергии, и можно ввести вений для квадратной решетки температуру. Обозначим число покоящихся частиц n0, частиц, имеющих единичную скорость n1, частиц со скоростью 2 n2. Тогда в каждом узле плотность = n0 + n1 + n2, полная энергия E = P + u2 /2 = = ni vi /2 = n1 /2 + n2 (P давление), температура T = P/. Так можно моделировать течения с переменной температурой, теплопередачу и выделение энергии. Несколько примеров столкновений, в том числе с выделением энергии, приведены на рис. 2.8.

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

2.4.2. Решеточное уравнение Больцмана В методе LGA существует значительный статистический шум, возникающий из-за случайности при столкновениях и из-за случайного начального распределения. Этот шум может быть полезен при моделировании процессов, в которых существенны флуктуации. Однако, для получения значений физических величин (плотности, скорости, давления, температуры) от него приходится избавляться путем усреднения по некоторой области пространства и/или за несколько шагов по времени. В случае детерминированных течений возможно также произвести несколько расчетов с одинаковыми (в среднем) начальными данными, а затем провести усреднение по всем экземплярам системы. Если число экземпляров стремится к бесконечности, то средние значения заселенностей называются одночастичными функциями распределения, которые будем обозначать fk. При отсутствии корреляций между одночастичными функциями распределения (гипотеза Л. Больцмана о молекулярном хаосе) эволюция системы описывается уравнением Больцмана. Так как у нас все значения заданы на решетке с небольшим числом возможных значений скорости и в дискретные моменты времени, имеем решеточное уравнение Больцмана, или LBE (Lattice Boltzmann Equation) [26]:

За шаг по времени все частицы должны перелететь в соседние узлы. Для этого возможные скорости частиц ck должны удовлетворять условию ck t = ek, где ek векторы, соединяющие узел с соседними. Обычно принимается t = 1. Левая часть уравнения соответствует перелету частиц из узла в узел (распространение). Последнее слагаемое отвечает за изменение заселенностей в столкновениях.

Макроскопические параметры плотность и скорость вещества в узле вычисляются по функциям распределения Очень простой и вместе с тем физически обоснованный вид столкновительного члена, предложенный в [27], описывает просто релаксацию системы к равновесному состоянию. Эта идея впервые была использована для метода LBE в [28] в виде:

Равновесные функции распределения fk зависят от плотности и скорости вещества в узле, чтобы выполнялись законы сохранения массы и импульса в столкновениях, то есть Кинетическая температура LBE частиц в энергетических единицах = kT /m задается уравнением и в простых моделях считается постоянной. Часто масса LBE частиц принимается за единицу, m = 1.

Такая модель хорошо описывает течения вязкой жидкости в пределе малых скоростей (число Маха M = u/cs 1, где cs = изотермическая скорость звука). Время релаксации определяет кинематическую вязкость = ( 1/2)c2 t. На твердых границах можно просто развораs чивать скорости прилетевших частиц, так моделируются непроницаемые стенки без проскальзывания.

* Приложение. Явный вид функций f eq. Обычно равновесные функции распределения выбираются в максвелловском виде:

В изотермических моделях достаточно разложить экспоненту в ряд с точностью до членов порядка u2, используя приближенную формулу В результате получаем:

Коэффициенты wk exp(ck /2) зависят и возможные векторы скорости только от модуля |ck |.

Необходимо подставить формулы для fk в уравнения (2.7, 2.8) и приравнять коэффициенты при всех степенях u вплоть до второй. Для одномерной модели c0 = 0, c1 = h/t, c1 = h/t получаем = 1/3(h/t)2, w0 = 2/3, w±1 = 1/6. В итоге, равновесные функции распределения имеют вид Здесь u = ut/h безразмерная скорость вещества.

В двумерной модели на квадратной сетке с 9 направлениями (рис. 2.9) c0 = 0); ck = h/t(cos(k/2), sin(k/2)) для k = 1... 4;

ck = 2h/t(cos((k + 1/2)/2), sin((k + 1/2)/2)) для k = 5... 8.

Здесь добавляется еще уравнение w(0)/w(1) = w(1)/w( 2). Получаем = 1/3(h/t)2, w0 = 4/9, w14 = 1/9, w58 = 1/36. Равновесные функции распределения в этом случае имеют вид где a = (t/h)2 / = 3, b = (t/h)4 /22 = 9/2, d = (t/h)2 /2 = 3/2.

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

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

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

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

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

2.4.4. * Модель LBE с внешними силами и фазовыми Метод LBE позволяет достаточно просто моделировать течения с фазовыми переходами жидкость – пар, причем возможно задать произвольное уравнение состояния. Кроме того, небольшое изменение уравнений позволяет рассчитывать течения, в которых на вещество действуют внешние силы.

Действие внешних сил Рассмотрим вначале, как моделировать силы, действующие на вещество.

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

При учете действия объемных сил решеточное уравнение Больцмана принимает вид:

То есть, после действия оператора столкновений (в котором используется скорость u) необходимо учесть изменение функций распределения fk под действием сил. Эта добавка равна разнице равновесных функций распределения при одной и той же плотности, но с разными скоростями [31] Можно изменить порядок сначала учесть действие сил, то есть вычислить промежуточные значения функций распределения а затем применить оператор столкновений Конечный результат будет тот же.

Часто нам необходимо знать скорость вещества (например, для расчета переноса примесей). Видно, что на каждом шаге по времени в каждом узле существуют два значения скорости до и после действия сил. Исследование уравнений LBE показывает, что в случае действия сил физическая скорость вещества равна их среднему арифметическому, то есть Фазовые переходы Достаточно простой способ моделирования фазовых переходов жидкость – пар был впервые предложен в статье [32]. Между частицами, находящимися в соседних узлах, задается сила взаимодействия:

Значения коэффициентов Gk > 0 соответствуют притяжению между соседними узлами, что необходимо для сосуществования жидкой фазы и паровой фазы. В обратном случае при Gk < 0 отталкивание. Здесь Gk выбираются таким образом, чтобы сила была достаточно изотропной (чтобы, например, капли получались круглыми). При использовании модели LBE на квадратной сетке сила взаимодействия между узлами, расположенными по диагонали на расстоянии 2, должна быть в 4 раза меньше, чем между ближайшими соседями, то есть G14 = G0 > 0, а G58 = G0 /4. Эффективная плотность () может выбираться достаточно произвольно.

Введение такого взаимодействия приводит к уравнению состояния среды в виде Коэффициент зависит от используемой решетки. Так, для одномерной модели с тремя значениями скорости частиц = 1, а в двумерной модели с 9 направлениями = 3/2. Напомним, что = 1/3(h/t)2 для рассматриваемых одномерной и двумерной моделей.

Теперь понятно, как ввести произвольное уравнение состояния P P (, T ) [33,34]. Введем новую эффективную плотность (, T ) = = G0 (, T ), которую выразим из уравнения (2.10):

Сила, действующая на вещество в узле, вычисляется в одномерном случае по уравнению, аналогичному (2.9) 2.4. Решеточные газы, решеточное уравнение Больцмана а в двумерном случае по уравнению где (x) = ((x)).

Здесь опять Gk = G0 для основных направлений и Gk = G0 /4 для диагональных.

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

F (x) = [A(x + 1) + (1 2A)(x) + A(x 1)] · [(x + 1) (x 1)] в одномерном случае и в двумерном. Коэффициент A зависит от используемого уравнения состояния.

Одним из достаточно простых уравнений состояния, позволяющих описать как жидкую, так и газообразную фазы, является уравнение Вандер-Ваальса (2.5). Удобно переписать его в безразмерных переменных Для согласования размерностей в методе LBE в правой части уравнения (2.11) необходимо давление P (, T ) домножить на коэффициент kef f = Pc /c · (t/h)2, который для инертных газов равен примерно k 0,01 (если взять h = 1 мм, t = 1 мкс такие единицы часто используются в физике быстропротекающих процессов). Для уравнения состояния (2.14) наилучшее значение параметра A в формуле (2.13) равно A = 0,152.

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

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

где индекс s означает номер вещества (1 или 2). Плотности и скорости каждого вещества в узле вычисляются по обычным формулам Общая плотность в узле равна сумме плотностей компонентов а выражение для общей скорости получается более сложным Такая формула необходима для сохранения импульса после действия оператора столкновений.

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

2.4.6. Химические реакции в решеточных моделях В статье [35] рассмотрено применение метода LGA для моделирования системы реакция-диффузия. Такие системы обладают многими интересными свойствами. Например, могут происходить колебания концентраций во времени. Хороший пример такого поведения дает реакция Белоусова – Жаботинского окисление малоновой или лимонной кислоты 2.4. Решеточные газы, решеточное уравнение Больцмана бромат-ионом, катализируемое ионами церия (про эту красивую реакцию можно прочитать, например, в [36]). В других случаях возможна самоорганизация возникновение устойчивых пространственно-неоднородных структур. Уравнение для системы реакция-диффузия в покоящейся среде имеет вид:

слагаемое C({n}) отвечает за химическую реакцию и обычно выражается в виде суммы произведений вида ks1 s2... na1 na2...

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

Гораздо более изящно описывается реакция в методе LBE [37]. Оператор столкновений делится на две части: гидродинамическую, обычного вида, и химическую, соответствующую уравнению реакции. Распространение и гидродинамическая часть оператора столкновений отвечают за диффузию, химическая часть за реакцию. Химическая часть записывается в том же виде, что и C({n}), необходимо только принять меры к тому, чтобы реакция не приводила к изменению скорости вещества.

* Пример. Брюсселятор. Приведем пример системы химических реакций:

Эта система называется Брюсселятором [38] в честь города Брюсселя, в котором работает предложившая данную модель группа Нобелевского лауреата И. Пригожина. Концентрации веществ A и B поддерживаются постоянными, вещества D и E удаляются из системы. Обозначим безразмерные концентрации веществ A, B, X, Y буквами a, b, x, y, соответственно. Дифференциальные уравнения для концентраций веществ X и Y выглядят следующим образом:

Для однородного состояния диффузия отсутствует, так как 2 x = 2 y = = 0. В этом случае равновесие в системе (2.15) достигается при x0 = a, y0 = b/a. Однако при b > a2 + 1 это равновесное состояние становится неустойчивым, и в системе возникают колебания концентрации веществ X, Y периодически изменяются. При усилении неравновесности возникают также колебания в пространстве образовываются волны концентраций, периодически пробегающие по системе.

Другой пример системы, обнаруживающей интересное поведение, можно найти в статье [37].

Достаточно полный (хотя и несколько устаревший) обзор работ, посвященных методу LBE, можно найти в статьях [39, 40], а также в Интернете по адресу:

http://www.scholarpedia.org/article/Lattice_Boltzmann_Method.

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

2. Рассмотрите одномерную задачу распада разрыва. Для этого в некотором вертикальном слое x1 < x < x2 первоначально задается повышенное значение плотности. На границах этого слоя плотность вещества испытывает разрыв. В одномерном случае все величины не зависят от вертикальной координаты. Постройте графики плотности и горизонтальной скорости в зависимости от горизонтальной координаты. Для этого их надо усреднить по y, то есть по вертикальным рядам узлов.

3. Смоделируйте разлет газового облака (круглой области), плотность которого выше, чем у окружающей среды.

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