Решение задачи#

Чтобы решить задачу (19)-(23) можно воспользоваться методом Ньютона-Рафсона. Однако, для изохорно-изотермической задачи двухфазного равновесия показана сходимость метода простой итерации. Единственное, вблизи критической точки сходимость медленная. Поэтому, как правило, используют комбинацию метода простой итерации и метода Ньютона для лучшего быстродействия.

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

\[\mu^{id}_{i}(T) + RT \ln[y_{i} \varphi_{i}(\mathbf{y}) P] = \mu^{id}_{i}(T) + RT \ln[x_{i} \varphi_{i}(\mathbf{x}) P].\]

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

(24)#\[y_{i} \varphi_{i}(\mathbf{y}) = x_{i} \varphi_{i}(\mathbf{x})\]

для каждого компонента \(i\). Введём обозначение

\[K_{i} = \frac{\text{мол. доля компонента $i$ в паре}}{\text{мол. доля компонента $i$ в жидкости}} = \frac{y_{i}}{x_{i}}.\]

Величину \(K_{i}\) называют константой равновесия (и, реже, K-фактором).

Преобразуем равенство (24). Во-первых, воспользуемся константами равновесия \(K_{i}\), во-вторых, чтобы сформулировать задачу о неподвижной точке, представим составы фаз как функции \(K_{i}\)

(25)#\[K_{i} = \frac{\varphi_{i}(\mathbf{x}(\mathbf{K}))}{ \varphi_{i}(\mathbf{y}(\mathbf{K})) }.\]

Полученная система является системой на неподвижную точку вида \(\mathbf{K} = g(\mathbf{K})\).

Осталось определить зависимости \(x_{i}(\mathbf{K})\) и \(y_{i}(\mathbf{K})\). Воспользуемся уравнением (23), учтём в нём \(L = 1 - V\) (22), \(y_{i} = K_{i} x_{i}\), и получим

\[x_{i} (1 - V) + (K_{i} x_{i}) V = z_{i} \iff x_{i} (1 - V + K_{i} V) = z_{i}.\]

Откуда получаем выражения составов фаз через константы равновесия \(K_{i}\) и мольную долю паровой фазы \(V\)

(26)#\[\begin{split}\begin{align} x_{i} &= x_{i}(K_{i}, V) = \frac{z_{i}}{1 + (K_{i} - 1) V}, \\ y_{i} &= y_{i}(K_{i}, V) = K_{i} x_{i} = \frac{z_{i} K_{i}}{1 + (K_{i} - 1) V}. \end{align}\end{split}\]

Чтобы получить зависимость составов фаз \(\mathbf{x}\) и \(\mathbf{y}\) только от констант равновесия \(\mathbf{K}\), достаточно найти связь \(V(\mathbf{K})\). Для этого рассмотрим разность (20) и (21)

\[\sum_{i} x_{i} - \sum_{i} y_{i} = 0, \iff \sum_{i} (x_{i} - y_{i}) = 0.\]

Воспользуемся (26) и получим нелинейное уравнение на мольную долю паровой фазы

(27)#\[F(V) = \sum_{i} \frac{z_{i} (1 - K_{i})}{1 + (K_{i} - 1) V} = 0.\]

Это уравнение называется уравнением Рэкфорда-Райса (Rachford-Rice). Оно связывает (совокупный) состав смеси \(\mathbf{z}\), мольную долю паровой фазы \(V\) и константы равновесия \(K_{i}\). Другими словами, уравнение (неявно) определяет искомую связь \(V(\mathbf{K})\).

В итоге имеем систему на неподвижную точку (25), где \(\mathbf{x}(\mathbf{K})\) и \(\mathbf{y}(\mathbf{K})\) определяются по константам равновесия через уравнение Рэкфорда-Райса (27) и уравнения (26).

Алгоритм решения#

По существу, алгоритм решения задачи сводится к построению системы на неподвижную точку (25). Пользователь сообщает давление \(P\), температуру \(T\), состав смеси \(\mathbf{z}\) и начальное приближение задачи \(\mathbf{K}^{0}\).

Один шаг алгоритма решения представляет собой простую итерацию, в которой…

  1. Решаем уравнение Рэкфорда-Райса (27) для \(\mathbf{z}\) и \(\mathbf{K}^{0}\). Из решения получаем мольную долю паровой фазы \(V^{0}\).

  2. Определяем составы паровой \(\mathbf{y}^{0}\) и жидкой \(\mathbf{x}^{0}\) фаз из уравнений (26) для \(\mathbf{K}^{0}\) и \(V^{0}\).

  3. Рассчитываем летучести паровой \(\varphi_{i}(\mathbf{y}^{0})\) и жидкой фазы \(\varphi_{i}(\mathbf{x}^{0})\).

  4. Проверяем равенства (25) в некоторой норме.

    • Если равенства удовлетворены, то решение задачи найдено.

    • Если равенства не удовлетворены, то получаем новое приближение констант равновесия \(\mathbf{K}^{1}\) по правилу простой итерации

    \[ K^{1}_{i} = \frac{ \varphi_{i}(\mathbf{x}^{0}) }{ \varphi_{i}(\mathbf{y}^{0}) }\]

    и переходим на шаг 1.

Начальное приближение#

Начальное приближение для алгоритма построим по результатам проверки устойчивости однофазного состояния. Если было проведено несколько попыток проверки, выбираем ту, для которой значение tangent plane distance (3) самое отрицательное. Обозначим состав (в мольных долях) фазы-зародыша из этой проверки как \(\mathbf{a}\). Если проверка предполагала газовое однофазное состояние, то примем \(\mathbf{a}\) как состав капли конденсата, в этом случае

\[K^{0}_{i} = \frac{z_{i}}{a_{i}}.\]

Если же проверка предполагала жидкое однофазное состояние, то примем \(\mathbf{a}\) как состав пузырька пара, в этом случае

\[ K^{0}_{i} = \frac{a_{i}}{z_{i}}.\]

Возможны и другие начальные приближения, например, для нормальных алканов неплохо работает приближение на основе закона Рауля (т.е. используется предположение об идеальных законах смешения). В серии расчётов возможно использование близкого по условиям и посчитанного ранее решения, т.е. решение в (\(P\), \(T\)) используется как начальное приближение решения в (\(P + \delta P\), \(T + \delta T\)).