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

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

Подытожим, что, по существу, задача свелась к решению нелинейной системы (6) для вектора неизвестных \(\boldsymbol{X}\).

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

  • состав в мольных долях (а не количество вещества)

  • давление

  • «тепловую энегию» \(R T\) (температура встречается реже)

Для вычисления коэффициента летучести удобно пользоваться вторым (безобъёмным) вариантом формулы (12).

В качестве солвера нелинейной системы мы возьмём метод Ньютона. Метод требует вычисления Якобиана (матрицы Якоби) системы, об этом будет указано ниже.

После решения системы остаётся вычислить значение функции \(F\) (3) в найденной точке. Вывод о стабильности термодинамического состояния оставим на пользователя алгоритма, поскольку неравенство \(F < 0\) требует учёта численных ошибок.

Якобиан#

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

Автоматическое дифференциирование это численный способ подсчёта значения функции и её производной (градиента, Якобиана, Гессиана,..) в заданной точке. Один вид автоматического дифференциирования, дифференциирование вперёд, реализован в пакете ForwardDiff.jl.

Например, решается система \(\boldsymbol{F}(\boldsymbol{x}) = \boldsymbol{0}\) из двух уравнений

\[\begin{split}\boldsymbol{F}(\boldsymbol{x}) = \begin{bmatrix} x_1^2 + 3 x_1 x_2 \\ \sin x_2 + x_1 x_2 \end{bmatrix}\end{split}\]

её точный Якобиан имеет вид

\[\begin{split}\boldsymbol{J}(\boldsymbol{x}) = \begin{bmatrix} 2 x_1 + 3 x_2 & 3 x_1 \\ x_2 & \cos x_2 + x_1 \end{bmatrix}\end{split}\]

Вычисление по аналитической формуле и дифференциированием вперёд дают одинаковые результаты для точки \(\boldsymbol{p} = [1, 2]^\top\).

julia> using ForwardDiff

julia> F(x) = [
    x[1]^2 + 3 * x[1] * x[2],
    sin(x[2]) + x[1] * x[2],
]

julia> J(x) = [
           2 * x[1] + 3 * x[2]   3 * x[1]         ;
           x[2]                  cos(x[2]) + x[1] ;
       ]

julia> p = [1, 2]
2-element Vector{Int64}:
 1
 2

julia> F(p), J(p)
([7.0, 2.909297426825682], [8.0 3.0; 2.0 0.5838531634528576])

julia> ForwardDiff.jacobian(F, p)
2×2 Matrix{Float64}:
 8.0  3.0
 2.0  0.583853

julia> J(p)[end, end]
0.5838531634528576

julia> ForwardDiff.jacobian(F, p)[end, end]
0.5838531634528576

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

В свою очередь, платой является скорость работы.

И хотя вывести Якобиан для кубического уравнения не слишком сложно, с начала 2000-х годов набирают популярность уравнения состояния из семейства SAFT. Можете оценить объём работы по выводу производных: см. первые и вторые производные свободной энергии Гельмгольца в приложении статьи [GS01]. Вероятно, кто-нибудь однажды сделает возможным задавать уравнения в исходном коде символьно, символьно их дифференциировать и компилировать эффективную версию (аналитическую реализацию).[1]

Грубое и уточнённое решение#

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

Чтобы увеличить отказоустойчивость (робастность, англ. robust) алгоритма, будем искать решение в два этапа. На первом этапе найдём грубое решение системы методом простой итерации с исходным начальным приближением. Соответствующая система на неподвижную точку (сравните с (6)) имеет вид

(15)#\[\ln{X_i^{k+1}} = h_i - \ln{\varphi_i(\boldsymbol{x}(\boldsymbol{X}^k))},\quad i=1,\dots,n.\]

В этой системе вектор неизвестных – \([\ln{X_1}, \dots, \ln{X_n}]^\top\).

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

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

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

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

Для углеводородов найдены более подходящие начальные условия. Мы воспользуемся следующим приближением [Mic82b]

\[\frac{\text{доля комп. в газе}}{\text{доля комп. в жидкости}} \equiv K_i = \frac{P_{c,i}}{P} \exp\big[5.42 (1 - T_{c,i}/T)\big]\]

где \(K_i\) – константа равновесия, \(P_{c,i}\) и \(T_{c,i}\) критические давление и температура компонента \(i\).

Начальных условий два и тестировать нужно оба

  • \(X_i = K_i z_i\) – образование пузырька пара в жидкости;

  • \(X_i = z_i / K_i\) – образование капли конденсата в газе.

Здесь мы пользуемся \(X_i\) как составом в мольных долях (см. (7)).

Алгоритму стабильности также понадобится сообщить типы фаз: тип фазы в однофазном состоянии и тип фазы-зародыша в двухфазном состоянии. Поскольку начальных приближений два, то всего может получиться до \(2^3 = 8\) вариантов. В простом случае достаточно воспользоваться двумя вариантами

  • Однофазный газ, в двухфазном состоянии фаза-зародыш – жидкость (капля конденсата);

  • Однофазная жидкость, в двухфазном состоянии фаза-зародыш – газ (пузырёк пара).

Собираем всё вместе#

Входные данные алгоритма

  • «Уравнение состояния»: параметры всех компонентов, правила смешения;

  • Состав смеси (в мольных долях \(\boldsymbol{z}\) или количества веществ \(\boldsymbol{N}\));

  • Давление смеси \(P\);

  • «Температура» смеси \(R T\);

  • Начальное приближение решения системы (6) \(\boldsymbol{X}^0\);

  • Пометка о типе однофазного состояния (газ или жидкость);

  • Пометка о типе фазы-зародыша (газ или жидкость).

Выходные данные алгоритма

  • Флаг «было ли найдено решение», т.е. сошёлся ли алгоритм;

  • Состав в мольных долях фазы-зародыша \(\boldsymbol{x}^*\);

  • Значение функции \(F\) (3) при найденном составе \(\boldsymbol{x}^*\).

Как упоминалось ранее, ответственность об определении стабильности по значению \(F\) ложится на пользователя алгоритма.

Этапы алгоритма

  1. Поиск грубого решения исходной системы методом простой итерации (15);

  2. Формирование системы (6) и её Якобиана для использования в методе Ньютона;

  3. Запуск метода Ньютона для системы (6) при начальном приближении, найдённом на шаге 1;

  4. Подсчёт состава фазы-зародыша \(\boldsymbol{x}^*\) (см. (7)) и \(F(\boldsymbol{x}^*)\) из найденного решения \(\boldsymbol{X}^*\).

Если на шаге 3 заданная точность решения не была достигнута, в выходных данных результат помечают как не найденный. Тем не менее, если значение \(F(\boldsymbol{x}^*) \ll 0\), то однофазное состояние можно трактовать как нестабильное.

Содержание