Метод Монте-Карло в предсказании термодинамических свойств углеводородов

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

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

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

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

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

1. Определение свойств тяжелых углеводородов;

2. Определение тепловых свойств природных газов;

3. Получение фазовых диаграмм, фазовых плотностей одно-, двух и многокомпонентных систем;

4. Определение адсорбционных свойств;

5. Определение растворимости газов, например, в полимерных покрытиях труб.

Итак, цель данной работы — получение термодинамических характеристик смесей углеводородов на основе моделирования методом Монте-Карло.

Были поставлены следующие задачи:

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

2. Рассмотреть интересные с точки зрения нефтегазовой промышленности бинарные системы, в том числе смеси, включающие в себя полимеры.

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

Для достижения целей и поставленных задач был использован модуль GIBBS, встроенный в интерфейс MedeA©, в различных ансамблях.

1. Метод молекулярного моделирования

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

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

1.1 Статистическая механика

Энтропия и температура

Согласно квантовой механике система может быть найдена в одном из собственных стационарных состояний полного гамильтониана. Для каждого их таких состояний |i> имеем, где H_гамильтониан системы, — энергия состояния |i>. Для рассмотрения систем с количеством частиц порядка 10²³ и огромной степенью вырождения уровней энергии необходимо применение аппарата статистической механики. Обозначим ? число собственных состояний системы с энергией E, состоящей из N частиц и занимающей объем V. Согласно постулату о равновероятности микроскопических состояний: система с фиксированными N, V и E с равной вероятностью может быть найдена в любом собственном состоянии ?.

Рассмотрим систему с полной энергией E, которая состоит из двух слабовзаимодействующих подсистем, то есть подсистемы могут обмениваться энергией, но полную энергию системы можно рассматривать как сумму энергий подсистем. При этом существует множество способов, которыми можно распределить полную энергию системы между подсистемами так, что E₁+E₂=E. И при определенном выборе Е₁ общее число вырожденных состояний системы ?₁ ?₂. Следует отметить, что общее число состояний не является суммой, но зависит от числа состояний в индивидуальных системах. Из этого следует, что удобно было бы иметь аддитивную меру вырождения подсистем. Поэтому логично использовать логарифм:

(1)

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

(2)

или другими словами:

. ()

Введя понятие

()

можно записать в виде

. (5)

молекулярный полимер механика инверсивный

Понятно, что поместив всю энергию в систему 1, будет происходить переход энергии из подсистемы 1 в подсистему 2 до тех пор, пока не выполнится условие уравнения, то есть подсистемы не окажутся в тепловом равновесии. При этом ln? всей системы максимален. Из этого можно сделать вывод, что она связана неким образом с термодинамической энтропией S. И существует несколько способов, по которым эти понятия могут быть связаны. Рассмотрим случай, когда

()

где — константа Больцмана.

Из условия теплового равновесия между подсистемами следует, что в₁= в₂, что другими словами: два тела находятся в тепловом равновесии, если их температуры равны. Из этого следует, что в связана с абсолютной температурой. Термодинамическое определение температуры:

()

Используя это, имеем, что:

. (8)

Введя понятие температуры, рассмотрим систему A, находящуюся в тепловом равновесии с тепловым резервуаром B. Полная система замкнута, поэтому полная энергия фиксирована. Предположим, что система А находится в квантовом состоянии i с энергией E_i. Тогда энергия резервуара E_B=E-E_i, а вырождение ?_B. Вероятность найти систему А в состоянии i:

. (9)

Для вычисления разложим вблизи 0:

(10)

или, используя выражения и:

. (11)

Подставив это выражение в, имеем:

. (12)

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

(13)

где при этом определяем статистическую сумму Q. Сравнив это выражение с термодинамическим соотношением:

где F - свободная энергия Гельмгольца. F связана со статистической суммой Q:

. (14)

Связь между свободной энергией Гельмгольца и статистической суммой зачастую более удобна в использовании, чем связь между ln? и энтропией. И как следствие последнее уравнение является рабочим уравнением равновесной статистической механики.

Классическая статистическая механика

В предыдущем параграфе нами была рассмотрена квантово-механическая формулировка основных положений и определений статистической механики. Энтропия связана с плотностью состояний системы с энергией E, объёмом V и числом частиц N. А энергия Гельмгольца связана со статистической суммой Q. Рассмотрим некую наблюдаемую А. Поскольку мы знаем вероятность того, что система при температуре Т будет обнаружена в собственном состоянии с энергией Е_i, то можем определить среднее тепловое значение наблюдаемой:

()

где — среднее значение оператора А в квантовом состоянии i. Таким образом, для вычисления средних сперва нужно решить уравнение Шредингера для системы, а затем посчитать среднее от оператора А для всех квантовых состояний со статистическим весом, которым нельзя пренебречь. Но, к сожалению, не для любой произвольной многочастичной системы можно решить уравнение Шредингера. И даже если решить уравнение возможно, число квантовых состояний, которые влияют на среднее значение, будет настолько велико, что численное вычисление средних будет невозможно. Однако, уравнение может быть упрощено в классическом пределе. Для начала перепишем это уравнение в том виде, когда оно не будет зависеть от базисного набора. Используя то, что:

где H — гамильтониан системы, можно написать:

(16)

где Tr — след оператора. Поскольку след оператора не зависит от выбора базиса, можно вычислять средние, используя любой базис. Предпочтительно использовать простой базисный набор, например, набор собственных функций оператора координаты или импульса. Также используем тот факт, что гамильтониан состоит из двух частей, соответствующих кинетической и потенциальной энергиям. Оператор кинетической энергии — это квадратичная функция импульсов всех частиц. Как следствие, собственные состояния импульсов также собственные функции оператора кинетической энергии. Аналогично, оператор потенциальной энергии — функция от координат частиц. Матричные элементы оператора потенциальной энергии удобнее всего считать в базисе собственных функций координат. Однако, H=K+U не диагонален ни в одном из перечисленных базисов. Но, если произвести приблизительную замену exp на exp exp, то можно упростить:

exp exp=exp {-в},

где — коммутатор операторов кинетической и потенциальной энергий. Легко проверить, что данный коммутатор пропорционален постоянной Планка ћ и будет пренебрежимо мал в квазиклассическом приближении ћ>0. Таким образом, в пределе получаем:

Tr exp?Tr exp exp

Обозначив |r> и |k> соответственно собственный вектор оператора координат и импульса, последнее выражение можно переписать:

Tr exp =

Все матричные элементы могут быть посчитаны напрямую:

,

где — функция координат всех N частиц. Аналогично,

,

где p_i=ћk_i, и

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

(17)

где d_размерность системы. Фактор 1L, но величина дисперсии может быть существенно уменьшена выбором w таким, что fw постоянной, при этом в целом дисперсия пропадет. В противоположность, если w постоянна, как в случае обычного отбора Монте-Карло, то относительная погрешность I может оказаться очень большой. Поскольку подынтегральная функция в не равна нулю только для тех конфигураций, где гиббсовский фактор ненулевой, то рекомендуется проводить неравномерный отбор Монте-Карло конфигурационных пространств таких, что весовая функция w примерно пропорциональна гиббсовскому фактору. Но, к сожалению, метод существенной выборки, описанный выше, не может быть использован для вычисления многомерных интегралов по конфигурационным пространствам, как в уравнении. Причина в том, что неизвестно, как построить такое преобразование из в, которое позволило бы генерировать точки в конфигурационном пространстве с плотностью вероятности, пропорциональной гиббсовскому фактору. Необходимое условие для решения этой проблемы — это возможность посчитать аналитически статистическую сумму исследуемой системы. Если бы мы могли это сделать, то вряд ли была бы необходимость в компьютерном моделировании.

Метод Метрополиса

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

. ()

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

Обозначим конфигурационную часть статистической суммы Z:

. ()

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

Предположим, что каким-то образом можно генерировать точки в конфигурационном пространстве согласно распределению вероятности. Это означает, что в среднем число точек n_i генерированных в единице объема вокруг точки, равно, где L — общее число сгенерированных точек. Другими словами:

. ()

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

Рассмотрим схему Метраполиса. Пусть — вероятность перехода из конфигурации o в n; m — число точек конфигурации o. Мы хотим, чтобы в среднем m было пропорционально N. Матричные элементы должны удовлетворять очевидному условию: они не должны разрушать равновесного распределения, если таковое установилось. Это означает, что в равновесии среднее число принятых пробных перемещений из системы о должно быть равно числу принятых пробных перемещений из n в о, то есть:

. (31)

Обозначим — вероятность совершить пробное перемещение из о в n, где обычно называют основной матрицей цепи Маркова. Обозначим вероятность принятия пробного перемещения из о в n. Таким образом:

. ()

Если симметрична, т. е., то можно переписать:

. (33)

Откуда следует, что

. ()

Существует множество вариантов, удовлетворяющих этому условию. Согласно Метрополису:

. (35)

Суммируя, в схеме Метрополиса вероятность перехода из о в n:

(36)

.

Предположим, что было сгенерировано пробное перемещение с U>U. Согласно уравнению это пробное перемещение будет принято с вероятностью:

<1. (37)

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

Необходимо отметить другое условие, что должно удовлетворять условию эргодичности.

Алгоритм метода Монте-Карло

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

1. Случайным образом выбирается частица и подсчитывается ее энергия U().

2. Частица перемещается на небольшое расстояние в случайном направлении r'=r+Д. Это приводит к изменению потенциальной энергии системы на некоторую величину ДU. Энергия частицы теперь U().

3. Перемещение из в принимается с вероятностью

. ()

Существует несколько типов перемещений в методе МК [9]:

1) Трансляция молекулы;

2) Вращение молекулы вокруг случайно выбранной оси;

3) Изменение объема;

4) Перемещение молекулы из одной ячейки в другую;

5) Уничтожение существующей молекулы или введение новой молекулы;

6) Увеличение какой-либо части молекулы;

7) Флип, то есть вращение одного атома вокруг оси, образованной его непосредственными соседями;

8) Рептация, то есть удаление одного конца молекулы и рост другого конца;

9) Перестановка, то есть удаление одной молекулы и помещение другой молекулы на ее место;

10) Вращение части молекулы вокруг некоторого атома;

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

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

1.2 Статистический ансамбль

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

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

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

, (39)

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

. (40)

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

Моделирование с использованием канонического ансамбля происходит по следующей схеме:

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

2. Частица перемещается на случайное расстояние,

,

где — максимальное смещение. Новая конфигурация обозначается n и ее энергия U.

3. Перемещение принимается с вероятностью:

.

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

(41)

где — число Авогадро, N — полное число частиц системы.

Изотермо-изобарический ансамбль NPT, в котором вместо объема постоянным поддерживается давление. Этот ансамбль очень часто используется, поскольку большинство экспериментов проводятся при контролируемых давлении и температуре. Предположим, что рассматриваются N идентичных атомов. Статистическая сумма имеет вид:

. (42)

При этом предполагается, что система находится в кубической ячейке с диаметром L=V^1/3, а r_i=Ls_i.

Отбор может производиться согласно правилу:

Более подробное описание может быть найдено в.

Рассмотренные выше ансамбли имели постоянное число частиц, но иногда необходимо получить информацию о среднем числе частиц в зависимости от внешних условий, например, в исследовании адсорбции в пористых твердых телах. Для таких целей целесообразно использовать Большой Канонический ансамбль, в котором постоянными параметрами являются температура, объем и химический потенциал µ_j j_го типа молекул. Для этого ансамбля статистическая сумма имеет вид:

(43)

А соответствующая плотность вероятности:

.

Вероятность принятия перемещения:

.

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

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

.

Итак, каждое микросостояние представляет собой моделируемую ячейку, содержащую исследуемый пористый материал и конфигурацию адсорбируемых молекул, определяемую T и µ. При этих условиях БКА позволяет флуктуировать плотности и энергии, при этом происходит выборка микросостояний, и подсчитываются средние значения флуктуирующих параметров. В результате адсорбционная изотерма представляет собой зависимость плотности от химического потенциала при постоянной температуре. Генерация микросостояний основана на процессе Маркова, т. е. при любой данной молекулярной конфигурации следующая генерируется путем случайного включения, удаления или перемещения адсорбируемой молекулы. Если же молекулы не сферические, то все движения сопровождаются случайным вращением.

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

1.3 Энергия молекулярной системы

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

(44)

где — расстояние между силовыми центрами. Параметры потенциальной энергии Леннарда-Джонса описаны в таблице 2.

Модель потенциальной энергии, рассматривающая индивидуальные атомы как силовые центры Леннарда-Джонса, называется «All Atoms» или «все атомы». Если силовой центр находится в некотором атоме группы атомов, то такая модель называется «United Atoms» или «объединенные атомы». Третья модель «Anisotropic United Atoms» или «анизотропно объединенные атомы», в ней силовой центр может находиться в пространстве между некоторой группой атомов (Рисунок 3).

Электростатическая энергия определяется выражением:

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

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

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

Рисунок 3. Принцип внутримолекулярного потенциала Anisotropic United Atoms и потенциала United Atoms. В первом случае AUA силовой центр расположен вблизи геометрического центра системы, во втором случае AU силовой центр находится в ядре атома углерода

Рисунок 4. Параметры гибкой молекулы

Эти потенциальные энергии вычисляются согласно следующим формулам:

энергия сгиба, (45)

энергия кручения, (46)

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

Таблица 2. Параметры Ленннарда-Джонса для групп CH₄, CH₃, CH₂. -расстояние от углерода до силового центра

2. Результаты расчетов

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

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

2.1 Параметры моделирования

Расчеты были проведены на модуле GIBBS, встроенном в интерфейс MedeA©. Были использованы периодические граничные условия для моделируемой ячейки. Был использован потенциал Леннард-Джонса вместе с правилом смешивания Лоренц-Бертело. Радиус обрезания составлял половину длины ячейки. Заряды не учитывались, поскольку число молекул с электростатическими зарядами было очень мало, и их электростатической энергией взаимодействия было пренебреженно. Для тяжелых алканов также были использованы потенциалы изгиба и кручения. Для получения сходящихся результатов было использовано порядка 3*10⁷ шагов, причем для усреднения были использованы только последние 50% данных.

2.2 Определение плотностей

Природный газ Бавлинского месторождения состоит из 7 основных компонентов. Общее число молекул в моделируемой системе было принято равным 500. В таблице 3 представлены эти данные. Аналогичная таблица представлена и для другого Ишимбайского месторождения республики Башкортостан. Эти смеси отличаются не только процентным соотношением компонентов, но и композицией: в смеси Ишимбайского месторождения также содержится сероводород.

Для расчета плотностей было использовано однофазное моделирование в NPT ансамбле. Для углеводородов были выбраны силовые поля AUA, описанные в разделе 1.4: для СН₄ был использован СН₄-Moller потенциал, для алканов AUA потенциал.

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

На рисунке 6 представлены результаты расчетов плотностей смеси природного газа Бавлинского и Ишимбайского месторождения в зависимости от давления от 30 до 110 МПа при пластовой температуре принятой равной 463 К и 310 K. Кроме того, расчеты были выполнены и для чистого газа метана. Экспериментальные данные для метана из согласуются с результатами моделирования.

Таблица 3. Процентное содержание компонентов смеси газа Бавлинского месторождения республики Татарстан [36]

Таблица 4. Процентное содержание компонентов смеси газа Ишимбайского месторождения республики Башкортостан

Рисунок 5. Плотности каждого компонента природного газа Бавлинского месторождения при 273 К и 0.013 МПа Рисунок 6. Зависимость плотности метана и природных газов Бавлинского и Ишимбаевского месторождений

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

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

2.3 Термодинамические коэффициенты

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

:

(47)

(48)

коэффициент сжимаемости

,

где — энтальпия, — внутримолекулярная потенциальная энергия, — число Авогадро, — число молекул.

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

где энтальпия H=U^ext+U^int+K+PV отличается от, потому что включает кинетическую часть K.

, H^id=U^int+K+NkT - энтальпия идеального газа, являющаяся только функцией температуры, H^res= U^ext+ PV - NkT - остаточная энтальпия.

.

Молярные теплоемкости определяются следующим образом:

гдечисло Авогадро.

Остаточная молярная теплоемкость определяется из расчетов Монте-Карло

Для смеси можно определить из идеальных теплоемкостей чистых компонентов, где — молярная доля компонента. Идеальные теплоемкости могут быть определены экспериментально, или с помощью метода групповых вкладов, но не путем моделирования, поскольку Монте-Карло моделирование не учитывает кинетическую часть энергии. В таблице 5 представлены идеальные теплоемкости чистых компонентов и смеси при Т=463К.

Таблица 5. Идеальные теплоемкости компонентов смеси и посчитанная идеальная теплоемкость смеси

Согласно уравнениям были рассчитаны в изобарно-изотермическом ансамбле тепловое расширение, изотермическая сжимаемость, фактор сжимаемости и теплоемкости. Результаты представлены соответственно на рисунках 7−10 в диапазоне давлений до 110 МПа при Т=463 К и 310 K.

Рисунок 7. Рассчитанное тепловое расширение природных газов

Рисунок 8. Рассчитанный фактор сжимаемости природных газов

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

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

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

2.4 Коэффициент Джоуля-Томпсона и инверсное давление

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

(49)

гдетеплоемкость, — молярный объем, — изотермическая сжимаемость.

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

Для газов коэффициент Дж-Т меняет свой знак с положительного на отрицательный при давлениях 30−60 МПа. Но проблема состоит в том, что измерять коэффициент Дж-Т очень трудно при подобных экстремальных давлениях и температурах. Эффект Джоуля-Томпсона очень слаб по сравнению с тепловыми потерями при высоких давлениях. А уравнения состояния возможны [4], но их результаты не подтверждаются измерениями при высоких давлениях. Кроме того, некоторые конденсатные газы могу содержать тяжелые углеводороды, которые делают предсказания неточными.

Дифференциальный коэффициент Дж-Т был определен совместным использованием теплового расширения и теплоемкости. Для этого использовался флуктуационный метод, разработанный в Lagache и другими, а также потенциал межмолекулярного взаимодействия AUA [38], упомянутый в разделе 1.4.

Результаты моделирования дифференциального коэффициента Дж-Т представлены на рисунке 11, где показана зависимость от давления при температуре 463 К. Эти результаты можно подтвердить путем измерения инверсного давления, то есть давления, при котором коэффициент Дж-Т меняет свой знак. Рассчитанные и экспериментально определенные инверсные давления представлены в таблице 6. Для сравнения также представлены результаты расчетов месторождения San Joaquin Valley из.

Рисунок 11. Рассчитанный коэффициент Дж-Т, T=463 K. Данные о природном газе San Joaquin Valle из [6]

Таблица 5. Инверсное давление при Т=463 К

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

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

где H — энтальпия.

Таблица 6. Интегральный коэффициент Джоуля-Томпсона для Бавлинского месторождения для диапазона давлений 10−50 МПа

Итак, метод Монте-Карло позволяет проводить вычисления широкого PT диапазона для определения энтальпий H=H. В результате этого можно построить изоэнтальпийные линии и вычислить интегральный коэффициент Джоуля-Томпсона. Данный метод дает результаты, коррелирующие с данными, полученными флуктуационным методом.

2.5 Фазовые диаграммы

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

Зависимости давления насыщенных паров природного газа Бавлинского месторождения и метана и их фазовые диаграммы представлены соответственно на рисунках 13 и 14. Также были определены критические давление и температура для метана и для природного газа. Результаты моделирования для метана хорошо согласуются с экспериментальными данными из.

Что касается смеси природного газа, то его критические параметры, а точнее псевдокритические, представляют собой средневзвешенные критические константы отдельных компонентов смеси. И согласно этому правилу критические параметры исследуемого природного газа Tc = 284 K, Pc = 4.6 MПa. Столь существенное отклонение и негладкий характер кривых для смеси природного газа можно объяснить тем, что каждый компонент смеси меняет свою фазу при разных условиях, в результате чего в смеси присутствует разделение и сосуществование разных компонентов в разных фазовых состояниях. Другой фактор, объясняющий подобное поведение, это флуктуации плотности и энергии вблизи критических точек. В нашем случае ситуация осложнена еще и многокомпонентностью смеси.

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

Рисунок 13. Давление насыщенных паров природного газа Бавлинского месторождения

Рисунок 14. Фазовые диаграммы природного газа Бавлинского месторождения и метана Рисунок 15. Зависимость параметра К_i от температуры для каждого компонента

2.6 Фазовые диаграммы бинарных систем

С точки зрения индустрии основная причина использования моделирования при рассмотрении бинарных систем — это обойти подгонку характерных параметров взаимодействия при встрече с новой смесью. Говоря о фазовом равновесии равновесия смесей, чаще всего встречаются соединения: CO₂, CH₄, H₂S, СН₄О и другие. Фазовое поведение таких смесей важно знать при разработке месторождений, особенно при расширении проектов по восстановлению обратной закачки газа.

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

Фазовые диаграммы углеводородных смесей обычно имеют форму петли на диаграмме. Причем максимальное давление может оказаться выше, чем критические давления входящих в состав смеси компонентов. На рисунке 16 показана фазовая диаграмма смеси CO₂-n_декана. Для CO₂ молекулы были использованы потенциалы All Atoms, а для цепочки декана AUA потенциал. Было показано, что правило объединения оказывает влияние на расчеты фазового равновесия в случае систем с существенно разными диаметрами, как в случае диоксида углерода и легких алканов. В нашем случае было использовано правило Лоренца-Бертелота.

Рисунок 16. Фазовая диаграмма CO₂-n_декан при 477 K

Как видно из рисунка 16 использованные параметры позволяют получить фазовые диаграммы с удовлетворительной точностью.

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

2.7 Растворимость газов в полимерах

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

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

В настоящей работе для оценки растворимости воды и углеводорода в полиэтилене был использован ансамбль Гиббса при постоянном давлении. Модель полиэтилена состояла из 15 линейных цепочек со степенью полимеризации 100. Для полиэтилена был использован AUA потенциал [38], а для молекул диоксида углерода потенциал All Atoms. Было использовано два этапа расчетов для получения сходящихся результатов.

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

Таблица 7. Результаты моделирования n-C₁₀₀+вода и n-C₁₀₀+CO₂ при 433 К., HDPE — high density polyethylene)

Рисунок 17. Фазовая диаграмма CO₂-nC₁₀₀ при 433 K из моделирования в Гиббс ансамбле

2.8 Определение адсорбционных свойств

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

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

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

Были использованы следующие параметры моделирования: молекулы С₃H₈, C₄H₁₀, C₅H₁₂ рассматривались как гибкие, остальные молекулы — как жесткие. Структура цеолита, параметры Леннарда-Джонса и заряды были взяты из базы данных MedeA. Решетка цеолита рассматривалась как неподвижная, а примитивная ячейка была расширена на 222, симметрия ячейки P1. Было использовано правило Лоренц-Бартелота объединения различающихся групп и суммирование по Эвальду для расчета электростатической энергии. Для воды был использован потенциал SPS, для твердого тела O-silicalite, Si-no-lj, Na-zeolite.

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

Рисунок 18. Структура Фоязита Na56Y

Рисунок 19. Адсорбционные изотермы компонент смеси Бавлинского месторождения при 308 К.

Таблица 8. Фугитивности компонентов смеси Бавлинского месторождения при 300 K

Рисунок 20. Количество поглощенного вещества компонент смеси Бавлинского месторождения в Фоязите Na56Y при 308 К

На рисунке 19 представлены адсорбционные изотермы компонентов смеси Бавлинского месторождения при 308 К от 1 до 10 бар. Для С₃H₈, C₄H₁₀ и С₅H₁₂ величина адсорбции постоянна, а для остальных компонент растет с давлением. Сильнее всего поглощается этан.

При рассмотрении смеси были подсчитаны фугитивности всех компонентов при T=308 К. Фугитивность данного газа — такая функция давления и температуры, подстановка которой вместо давления в термодинамическом уравнения для идеального газа делает их справедливыми и для реального газа при рассматриваемых условиях. Используя полученные фугитивности, было посчитано количество адсорбированного вещества для каждого компонента.

Заключение

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

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

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