Метод молекулярной динамики для исследования движения частиц веществ

новесному подходу происходит снижение концентраций метанола и этанола при полном испарении до нуля.

t, ч

Рис. 3. Динамика изменения температуры на поверхности раздела фаз: Тг - температура по равновесной модели, Тп - температура на поверхности раздела фаз по неравновесной модели

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

ВЫВОДЫ.

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

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

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

УДК: 539.19:519.7

А.Л. Колобова, О.С. Кунина, Э. М. Кольцова

Российский химико-технологический университет им. Д.И. Менделеева, Москва, Россия

МЕТОД МОЛЕКУЛЯРНОЙ ДИНАМИКИ ДЛЯ ИССЛЕДОВАНИЯ ДВИЖЕНИЯ ЧАСТИЦ ВЕЩЕСТВ

The bases of molecular dynamics method with reference to solving numerical problems are stated in the given work. Various potentials of interaction are described; the parameters describing laws of interaction are determined; the basic assumptions which allow us to establish connection between microscopic parameters of modeling and macroscopic parameters of the modeled object are stated.

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

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

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

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

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

Перед тем как использование больших ЭВМ получило широкое распространение при расчете свойств атомных систем, необходимо было, чтобы потенциалы имели простую аналитическую форму. Наиболее популярными из таких стандартных потенциалов оказались экранированный кулоновский потенциал, потенциалы Борна-Майера, Леннарда-Джонса и Морзе[3].

Литературные данные показывают, что для многих оксидных систем ионная модель позволяет получить данные, которые неплохо согласуются с эксперименталь-

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

Применимость ионной теории для построения модели оксидных систем обуславливается значительной разностью электроотрицательностей кислорода (3.5 по Полингу) и большинства металлов (0.9 у Ка, 1.0 у Са, 1.5 у А1 и т.д.). Например, для СаО разность электроотрицательностей составляет ~2.5 и степень ионности связи Са-О близка к 76%. Для Б102 степень ионности составляет ~51%. Поэтому, несколько упрощая ситуацию, можно приписать частице кислорода формальный заряд - 2е (е - элементарный заряд), частицам Са, А1, - соответственно +2е, +3е, +4е. В случае ионной связи потенциал

взаимодействия двух частиц является парным и записывается обычно в форме потенциала Борна-Майера (1) или потенциала Полинга.

где г - расстояние между центрами ионов i и у, Zi, Z у - заряды /'-го и у-го ионов соответственно в единицах элементарного заряда е, В/у и Я/у -параметры, учитывающие отталкивание ионных оболочек, С/у и О/у -параметры, учитывающие диполь-

дипольное и диполь-квадриупольное взаимодействия соответственно

Впервые значения потенциалов Борна - Майера - Хаггинса (БМХ) и Полинга для близких к оксидам систем (галогенидов щелочных металлов) были определены в работах ¥.0.¥пш/, М.Р.Тоз/. Подгонку этих параметров проводили по таким известным экспериментальным характеристикам кристаллических соединений, как плотность, сжимаемость, энергия системы. По достигнутому уровню точности моделирования указанных свойств потенциал БМХ оказался лучше потенциала Полинга. В общем случае члены, учитывающие диполь-дипольные и диполь-квадрупольные взаимодействия, в выражении (1) играют большую роль. Например, для жидкого КС1 полная энергия при 1045К составляет 730.2 кДж-моль-1, а вклады от этих взаимодействий соответственно равны -26.5 и -4.1 кДж-моль- .В оксидах они еще менее существенны. Поэтому часто эти члены просто

Для кристаллических оксидов МО (М = М^, Са, Бг, Ва) параметры потенциалов Борна - Майера были рассчитаны по значениям параметров решетки и диэлектрическим постоянным в работе M.J.L.Sangster, А.М.&опекаш. Парные потенциалы для таких оксидов, а также для Б102, и02, ТЮ2, КЬ205, силикатов рассмотрены в работе С.ЯА.Сайо^, С.Ы.Егеешап. Параметры взаимодействия для пар М - О (М = Мп, Бе, Со, N1) оценены в исследованииЫ^опекаш, M.J.L.Sangster.

Несмотря на простоту, потенциалы (1) либо аналогичные потенциалы с добавками вкладов типа потенциала Морза можно использовать для описания кристаллических структур кварца, кристобалита, оксида бора, а также структуры жидких оксидов системы Ка20 - В203 и СаА12Б1208 рассмотрено в статье К. Ка^ашжа.

Если при моделировании не требуется рассчитывать энергию системы и связан-

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

ат

работ (например, в статье М.Атіпі, $.К.Міїга) использовали силовую функцию типа функции Полинга:

Выражение для потенциала Борна-Майера имеет вид:

иу (т)= + Ба ехР

опускают.

2

е

г + г

I ]

1 + sgn (г^]

Здесь sgn принимает значение + 1 или -1 в зависимости от знака аргумента (гг,), а

параметр п изменяется обычно в пределах от 8 до 10.

Авторы работы рассчитали потенциалы взаимодействия пар М - О (М = Са, В, А1, Б1, V) с помощью «магической» формулы Малликена, которая учитывает ионноковалентный характер связи. Полученные зависимости потенциалов от г аппроксимировали выражением Полинга:

и, (г )= гг

г

1 + ^п(гг]

\п

(п +1)-

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

* Оксид Zl ^2 г1, пм г2, пм

СаО 1.46 -1.46 102.0 133.0

В2О3 2.34 -1.56 10.0 140.0

А12О3 2.55 -1.70 42.5 140.0

Б1О2 2.30 -1.15 30.1 120.4

V2O5 3.65 -1.46 40.1 133.0

Нумерация атомов: 1 - М, 2 - О.

Авторами работы Л.И.Воронова, Б.М.Лепинских была проведена также оптимизация величин г и г по результатам квантово-химических расчетов тройных систем, включающих Са, В, А1, Б1, V и О. При этом получили

го =-1.46, = 2.92, Ху = 3.65 . В дальнейшем для расчета профилей ионно-

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

Потенциалы Борна-Майера использовались при моделировании структуры и свойств соединений ЫЫЬО3 и КпЬО3. Для всех пар ионов было принято, что

г + г] . Заряды и радиусы ионов подбирались так, чтобы полу-12.36

чить наилучшее согласие с наблюдаемыми фононными спектрами. Для кислорода в КпЬО3 пришлось выбрать два значения г, соответствующих двум различным позициям в решетке. Формальные заряды 2о сильно отличаются от обычно принятой величины -

2. Для системы Ы2О - В2О3 применялись парные потенциалы типа БМХ. Радиусы ионов равнялись 93 (Ы+), 74 (В3+) и 142 пм (О2-) а заряды 1,3 и -2 .

Описание межчастичного взаимодействия в ковалентном варианте проводится на основе трехчастичных потенциалов. Энергия системы записывается в виде:

В, =1822эВ

и

и

и о + С2 ЕЕ и

п

2

е

1

к

Здесь и0 - энергия идеального кристалла (идеальным считается кристалл, в котором длины всех связей и углы между связями равны равновесным значениям; равновесный угол определяется видом кристаллической решетки), С2, С3 - коэффициенты; функции и (г, Г] ) и у(т;, Г], Гк ) описывают изменение энергии из-за растяжения (сжатия) связей и

искажения валентных углов соответственно. Коэффициенты С2 и С3 подбираются таким образом, чтобы получить наилучшее согласие результатов моделирования и экспериментально установленных характеристик оксида (колебательные спектры и межъядерные расстояния). Наиболее простой вид функций и(г, Г] ) и у(г{ , Г], Гк ) был предложен Китингом:

u(г\,Г2 )= г - гГ - й02,

где Г и Г> - радиусы-векторы соседних ионов, й0 - равновесная длина связи,

АГі, Г] , Гк ) = [с^ в(ук) - СОБ во (і]к)]2,

где в]) - угол между связями І-] и ]-к, а во (і]к) - соответствующий равновесный угол. Суммирование в выражении (2) проводится по ближайшим соседям центральной частицы.

Парный потенциал Ми :

ab

и

,n rm

где 10 < n < 12 и 4 < m < 6;

Особенно часто трехчастичный потенциал применяли для кремния. Если в ионном варианте выбор начального состояния при моделировании мало существен (это может быть кристаллическая решетка или случайная система ионов), то в «ковалентном варианте» приходится специально «готовить» начальное состояние, в котором координационные числа фиксированы (например, 4 у Si и 2 у О). Такая подготовка начального состояния требуется, как правило, и при моделировании «рыхлых» систем обратным методом Монте - Карло.

По экспериментальным данным для SiO2 были подобраны коэффициенты для низко- и высокотемпературных кристобалита и кварца. Для пар Si-Si были использованы данные для чистого кремния и потенциал u(r) выбран в форме потенциала Леннарда-Джонса. В итоге удалось получить модель чистого кремния, устойчивого в решетке алмаза. Для пар О-О были использованы свойства молекул О2 и Оз (потенциал Леннарда-Джонса). При этом было получено соответствующее литературным данным значение энергии диссоциации кислорода, но стабильной конфигурацией молекулы О3 оказалась линейная, а не треугольная. Для пар Si-О также использовался потенциал Леннарда-Джонса. Остальные параметры модели подгоняли, добиваясь согласия между результатами расчетов и экспериментальными данными для кристаллического кремнезема. Рассчитанные ab initio потенциалы были использованы для описания параметров решеток и температур плавления оксидов MgO, Mg2SiO4 и MgSiO3.

В ряде исследований применялись комбинированные потенциалы для оксидов, включающие как ковалентные, так и ионные вклады. Моделировали стекла системы Na2O-SiO2, используя комбинацию потенциала Китинга для троек атомов, включающих Si и О, и потенциала Леннарда-Джонса для пар Na-Si и Na-О. Аналогичный прием применен для моделирования стекол типа Na2O-B2O3. Для троек О-В-О и О-В-В выбрали трехчастичный потенциал, а для пар, включающих Na, - парный потенциал Борна-Майера (1) [4].

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

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

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

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

Список литературы

1. Билер Роль машинных экспериментов в исследовании материалов// Машинное моделирование при исследовании материалов. - М.:Мир. - 1974. - с.31-250.

2. Псахье С.Г., Негрускул С.И. Зольников К.П. Дискретные компьютерные модели конденсированных сред с внутренней структурой // Физическая мезомеханика и компьютерное конструирование материалов: в 2 т. - Новосибирск: Наука. Сиб. изд. фирма РАН, 1995. -т.2. - с.77 - 99.

3. Каширин В.Б., Козлов Э.В. Влияние потенциала взаимодействия на структуру и свойства моделируемых аморфных структур // ФММ. -1993. - т.76, вып.1. - с.19-27.

4. Белащенко Д.К. Компьютерное моделирование некристаллических веществ методом молекулярной динамики.//Соровский Образовательный Журнал. 2001. Т. 7, №8. С. 44 -50.

УДК 681.39:663.1

Л.С. Гордеев, И.И. Дубровский, В.Л. Лукьянов, А.М. Шеховцова Российский химико-технологический университет им. Д. И. Менделеева, Москва, Россия

МИКРОПРОЦЕССОРНАЯ СИСТЕМА УПРАВЛЕНИЯ ПРОЦЕССОМ ФЕРМЕНТАЦИИ

In this article the main principles of construction of microprocessors management systems are stated. They were adopted for the biotechnology processes. The variant of microprocessors control systems on the base of programmed controllers for the device for fermentation is shown. Controllers make discrete signals, which control the dose of resources and materials during loading device, and carries out the subsequent management of process of a fermentation.

Изложены основные принципы построения микропроцессорных систем управления (МПСУ), адаптированные к процессам биотехнологии. Показан вариант создания МПСУ на базе программируемых контроллеров для ферментера, осуществляющие формирование дискретных сигналов, управляющих дозированием сырья и материалов при загрузке аппарата и осуществление последующее управление процессом ферментации.

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

Применяя современные принципы построения микропроцессорных систем управления (МПСУ), можно определить следующую структуру данной системы:

1) уровень датчиков и исполнительных механизмов;

2) уровень контроллеров и модулей ввода-вывода;