Оптимизационный алгоритм решения задачи о фазовом равновесии

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

В общем случае условия фазового равновесия должны быть согласованы с уравнениями состояния фаз многокомпонентной смеси [3,4]. В данной работе рассматривается приближение, когда, коэффициенты летучести не зависят от составов фаз (или этой зависимостью можно пренебречь). В этом случае условия фазового равновесия описываются в терминах коэффициентов распределения, зависящих только от давления и температуры [5]. Существует точка, зрения, что все компоненты в той или иной мере содержатся во всех фазах [6]. В данной работе, однако, считается, что часть компонентов не содержится в некоторых фазах. В этом смысле рассматривается ограниченная задача, о фазовом равновесии. Компоненты смеси разделяются на. два. типа: активные — те, которые участвуют в установлении фазового равновесия, и инертные, роль которых ограничивается участием в балансе энергии. Принимается, что активные компоненты могут образовать три фазы: жидкую (нефтяную), газовую и водную. Жидкая фаза, содержит компоненты, не присутствующие в других фазах (тяжелые нефти), и компоненты, присутствующие в газовой фазе (легкие нефти). Газовая фаза, наряду с легкими нефтями содержит нерастворимые газы (содержащиеся только в газовой фазе) и пары H2O. Водная фаза состоит только из одного компонента H2O. Количество фаз, на которое расслаивается смесь, заранее неизвестно и определяется термобарическими условиями и составом смеси. Инертные компоненты образуют скелет — набор неподвижных фаз.

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

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

Задача минимизации потенциала Гиббса (при фиксированных температуре, давлении и составе смеси) ставится в пространстве независимых степеней свободы: молярных долей (относительно полного количества молей) H2O в водной фазе и легких нефтей в жидкой фазе. Указанные величины однозначно описывают состояние системы (молярные доли и компонентный состав фаз). Условия минимизации потенциала Гиббса в такой постановке дают критерии фазового равновесия для сосуществующих фаз, а балансы количества молей компонентов (уже учтенные в аналитическом виде потенциала) приводят к уравнениям Речфорда-Райса для молярных долей фаз. Множество, на котором происходит минимизация ограничено и описывается набором естественных условий неотрицательности молярных долей компонентов.

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

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

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

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

Предлагаемый метод расчета фазового равновесия состоит в следующем:

• 1)    выбор моделвного потенциала Гиббса;

• 2)    модификация потенциала в соответствии с методом логарифмических барверов;

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

• 4)    добавление уравнения энтальпии, как неявного уравнения на температуру;

• 5)    совместное решение уравнений пунктов 3 и 4.

2.    Постановка задачи и допущения

В результате расчета фазового равновесия в и. 5 предлагаемого алгоритма активные компоненты смеси всегда образуют трехфазную систему. Если система фактически (т.е. с исходным модельным потенциалом Гиббса) находится в одно- или двухфазном состоянии, то молярные доли «паразитных» фаз не должны превышать некоторого порога, регулируемого высотой внутреннего барьера.

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

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

Рассматривается многокомпонентная система, которая может расслаиваться на несколько фаз. Часть из компонентов (инертные) не участвует в фазовых превращениях и формирует инертные фазы. Другие компоненты (активные) могут потенциально образовывать три подвижные фазы: жидкую (нефтяную) фазу (индекс L), газовую фазу (индекс G) и водную фазу (индекс ГТ).

Предполагается следующее распределение компонентов по фазам:

• •    Один компонент H2O — может присутствовать в водной и газовой фазах.

• •    Легкие нефти — могут присутствовать в жидкой и газовых фазах (например, C5H12 в примере, рассмотренном ниже).

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

• •    Нерастворимые газы — могут присутствовать только в газовой фазе.

Вводятся обозначения (все молярные доли рассматриваются по отношению к количеству молей всей смеси активных и инертных компонентов):

• •    с,- молярные доли H2O (г = 0), легких (г Е L) и тяжелых (г Е H) нефтей, нерастворимых газов (г Е G) II инертных компонентов (г Е S) в смеси:

• •    Х і ,у і — концентрации легких нефтей в жидкой и газовой фазах, уо — концентрация H2O в газовой фазе: концентрация H2O в водной фазе равна. 1:

• •    А, ш,7 — молярные доли жидкой, водной и газовой фаз в смеси соответственно. При этом А + ш + у = 1, т.к. учитываются также инертные фазы;

• •    в водной фазе находится только H2O, а значит, ш является молярной долей H2O в водной фазе;

• •    Аі = Ажі, г Е L — молярные доли (относительно полного состава смеси) легких нефтей в жидкой фазе.

В табл. 1 приведен пример распределения компонентов по фазам. Компоненты C5H12 и CioH22 считаются легкими нефтями, компонент C15H32 — тяжелой нефтью. В смеси присутствуют три нерастворимых газа (N2, CO2 и O2), а также два инертных компонента-фазы (Sand 11 Gran).

Таблица!

Пример распределения компонентов по фазам

\	L (нефть)	G (газ)	W (вода)	S (песчаник)	Gr (гранит)	∑︀
H2O		со — ш	ш			с0
C5H12	Х1	с1 — Х1				сі
C10H22	Х2	С2 — X2				с2
C15H32	сз					сз
N2		с4				с4
CO2		с5				с5
O2		сб				сб
Sand				с7		с7
Gran					с8	с8
∑︀	X	7	ш	с7	с8	1

Требуется найти количество фаз, на которое расслоится смесв, молярные доли этих фаз и их состав, при известных молярных долях компонентов сі и констант фазового равновесия К(р,Т),г Е L. Здесь неизвестными величинами являются молярная доля водной фазы ш и молярные доли легких нефтей в жидкой фазе X i — в сего ||L| + 1 неизвестная.

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

Принимается, что молярный потенциал рассматриваемой системы имеет вид

9 = 9S + X У^ Xi(m,L + НТlnж) + ieL,H

+ 7 У2 ^-G + НТlnУІ) + 7№(№,G + НТlnУо) + ^o,w, ieL,G где 9S ~ часть потенциала, обусловленная скелетом; р^(р,Т), рі,и(р,Т ),^o,w (р,Т) — химические потенциалы чистых компонентов в жидкой, газовой и водной фазах; Н — газовая постоянная.

После исключения из (1) несущественных для дальнейшего слагаемых (зависящих только от фиксированных давления р, температуры Т и состава с^ модельный потенциал Гиббса принимает вид у = шко +Е Xi^+EL(Xi) — L(x)+ ieL         ieL

+ L(cо — ш) +     L(c i — Xi) — L(s — X — ш).

ieL где Ki

X = Ex

^ i,L — ^ i,G НТ

= In Ki (р,Т) — логарифмы констант равновесия: L ( x ) = ж іпж:

i'. а молярная лол>i подвижных фаз s = 7 + X + ш = const. В общем случае

ieL,H s = 1, так как задача сформулирована для всех компонентов, а не отдельно для подвиж

ных фаз.

Минимум функции у ищется при естественных ограничениях:

{ 0 6 Ш 6 Со,

0 6 X i 6 c_i ,     г G L.

Если минимум потенциала (2) достигается строго внутри области, заданной ограничениями (3), то выполняются условия фазового равновесия:

Эр ЭХ ч Эр Эш

Хі/X      , Кі ц г G L,

Ki + ln 7----—— = In = 0,

(C i - X i )h

1        1 Ко

Ko + In  -----—- = In — = 0.

(co - w)/7yo

3.    Энтальпия

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

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

• •    зависимость фазовых составов от температуры может содержать разрывы первого рода.

4.    Метод решения

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

В рассматриваемой задаче молярная энтальпия принимается в виде

К = £ Ci hi (р, Т) - шАһо(р, Т) - £ Хі Ahi (р, Т), (4) i iGL где Ahi(р,Т) — удельная молярная энтальпия парообразования (положительная величина), a hi(р,Т) — молярные энтальпии компонентов в «основной» фазе. Для H2O и легких нефтей, которые участвуют в фазовых переходах, «основной» фазой является газовая G.

С другой стороны, энтальпия может являться не вспомогательным, а определяющим параметром. Условия фазового равновесия системы (минимизация функционала (2) по неизвестным { ш,Хі }) и соотношение (4) совместно рассматриваются как система нелинейных уравнений для определения температуры и соответствующего этой температуре фазового равновесия. Эту систему можно трактовать следующим образом: требуется подобрать такую температуру смеси, чтобы после установления фазового равновесия ее молярная энтальпия равнялась заданной величине.

Отметим, что условие (4) не является ограничением в задаче (2). Соотношение (4) — это уравнение на Т, а блок фазового равновесия — это неявная функция фазовых составов от температуры.

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

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

В соответствии с этими рассуждениями произведем коррекцию функционала р в (2) следующим образом:

ф = р — е

ln Ш + ln(co

^{— ш) +} ^Щ ( ^{ln A} i ^{+ ln(c}i ^{— A}i)) iCL

Задача минимизации функционала р заменяется задачей минимизации функционала ф. Локальный минимум функционала ф находится строго внутри области вблизи (при достаточно малых е и еф искомого минимума функционала р. Исходную задачу можно переформулировать в виде: найти минимум функции ф, заданного в (5), по переменным {ш,АД, при условии Һ ( ш,А і ,Т ) = ту, где у — заданное значение молярной энтальпии. Отметим, что последнее условие является не ограничением в задаче минимизации функционала ф (модельный потенциал Гиббса ф не минимизируется по температуре T), а дополнительным уравнением для определения температуры смеси. Таким образом получаем систему нелинейных уравнений:

0 = F ( x ) = <

РЩх) Эш Э-ф (х) ЭХ і

л(х) —

ту

где х = ( ш,А і ,Т ).

Система (6) далее решается модифицированным методом Ньютона:

A x k = — ( ^ F d x k² Г F < X ‘),

X k +1 = X k + « k ^x k , « k = min(1,3«k,max).

В методе (7) используется ограничитель шага « к 6 1 Для того, чтобы не выйти за границу области определения функции. Когда « к = 1, шаг происходит согласно обычному методу Ньютона. «k,_max > 0 — это максимальное число, такое что Xk + «k,_maxAxk еще принадлежит области ш Е [0, co], Ai Е [0,сф а также выполняется условие T Е [Tmin, T_max]. Множитель ф <  1 (типично ф = 0.9 ^ 0.99) используется, чтобы новое приближение было строго внутри области.

Решение задачи можно аналитически продифференцировать по параметрам:

Ғ(х(а), а) = 0, дF(x, а) дx   дF(x, а)

Эх  да    да дх д а

^д Ғ ( х , а ) ^ ¹ д Ғ ( х , а )

5.    Результаты

Для верификации предложенного алгоритма рассматривалась задача о нахождении фазового равновесия смеси, включающей H2O (молярная доля 0.5) и две легкие нефти: C5H12 и C10H22 (молярные доли 0.25). При этом задача решалась в диапазоне температур 530 660 К при давлении 107 Па. На рис. 1 показаны результаты расчета.

Как видно из графика, большую часть диапазона температур занимают двухфазные состояния (при T < 550 — водная и жид кая фазы, а при T > 560 — газовая и жидкая). Такие двухфазные состояния равновесия, как отмечалось выше, соответствуют граничным точкам области минимизации потенциала Гиббса и недостижимы в предложенном алгоритме. Однако визуально построенный график идентичен такому же графику, полученному из «стандартного алгоритма». Здесь под «стандартным» понимается рассмотрение случаев фазового состава смеси в состоянии равновесия с точным выполнением соотношений Кі(р,Т) = — для компонент, которые находятся в нескольких сосуществующих фазах.

Наибольшее отличие результата, полученного предложенным методом, от «стандартного» наблюдается в окрестности перехода от двухфазного к трехфазному состоянию (рис. 2). Однако и в этом случае разница невелика — она составляет не более 0.5%.

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

6.    Заключение

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

Исследование выполнено за счет гранта Российского научного фонда (проект № 15-11-00015) в лаборатории флюидодинамики и сейсмоакустики МФТИ. Авторы выражают благодарность А. А. Тараканову за предоставление «стандартной» версии алгоритма для сравнения.