Расчет фазовых диаграмм с TDLIB

http://evgenii.rudnyi.ru/doc/misc/tdlib.html

Файлы доступны в архиве http://evgenii.rudnyi.ru/soft/tdlib00+.tar.gz.

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

При компиляции библиотеки получается консольное приложение assess.exe, которое как раз может решать обратную задачу: определение энергии Гиббса из доступных экспериментальных данных. Однако, поскольку решение прямой задачи (расчет равновесного состава) является частью обратной задачи, то данное приложение можно использовать только для прямой задачи. Решение прямой задачи будет коротко рассмотрено в настоящем документе.

Должен отметить, что при расчете фазовой диаграммы наиболее оптимально использовать минимизацию энергии Гиббса (см. уравнение 9.1 в http://evgenii.rudnyi.ru/doc/ru/teaching/313/313group.pdf), например см. ПО для расчета равновесного состава в трехкомпонентной системе Ba-Cu-Y

http://evgenii.rudnyi.ru/soft/bacuy_eq.tar.gz

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

Соответственно TDLIB предлагает только возможность расчета равновесного состава при заданном составе фаз, то есть решение системы нелинейных уравнений 2.26 (см. документацию к библиотеке). Рассмотрим как это может быть использовано при расчете фазовых диаграмм на примере фазовой диаграммы Bi-Se из работы http://evgenii.rudnyi.ru/doc/ru/papers/01nm_bise.pdf.

Соответствующие файлы находятся в директории ex/bise. Ниже предполагается, что это текущая директория и программы assess.exe и gnuplot.exe (wgnuplot.exe) находятся на пути (PATH), заданном ОС.

Файл ex/bise/sys.mod задает мольные энергии Гиббса всех фаз, которые принимаются в расчет. Это исходная информация для расчета фазовых диаграмм в прямой задаче. Файл ex/bise/alg.mod описывает возможные фазовые равновесия в системе. Файл ex/bise/pd.out.mod описывает постороение фазовой диаграммы. В нем рассчитываются фазовые равновесия при заданном составе фаз. Здесь получается, что уже необходимо представлять себе возможную топологию диаграммы состояния и требуется задать расчет всех линий в диаграмме.

Команда

$ assess sys alg pd.out -o t

считывает энергии Гиббса, информацию о фазовых равновесиях и затем проводит расчет фазовой диаграммы, которая будет записана в файл t.pd в формате gnuplot. Команда

$ gnuplot t.pd -

построит рассчитанную фазовую диаграмму. Более подробное описание объектов в файлах ввода можно найти в документации к TDLIB (см. раздел 3.4.7. Object PhaseEquilibrium (tdlib/ex/bise), а также другие разделы).

Следующее

Введение в TDLIB: чистые вещества
Введение в TDLIB: регулярный раствор
Введение в TDLIB: объект фазовые равновесия
Уравнивание химических реакций в TDLIB
Программа VCS из книги Smith&Missen


Comments

9 комментариев to “Расчет фазовых диаграмм с TDLIB”

Comments are now closed
  1. Уважаемый Евгений,
    Здравствуйте!

    Прежде чем «окунуться» в расчет фазовых диаграмм с TDLIB
    в качестве вводного курса прочитал ваши лекции 313group.pdf (жаль, что мне не повезло учиться на хим. ф-те МГУ). Интересно написано.

    Как я понимаю, это первый вариант и, к сожалению, не все вопросы остались раскрытыми полностью (интересно было бы посмотреть литературные ссылки, по крайней мере [2, 5]). В этой связи появилось несколько вопросов:

    1.
    На с.14 для смеси компонентов H2, N2, NH3, N2H4
    «…протекает две химических реакций.»
    3 H2 + 1 N2 == 2 NH3
    2 H2 + 1 N2 == 1 N2H4

    Интересно почему рассматриваются только эти две кандидатуры?

    По моим расчётам,
    могут протекать как минимум 5 реакций, не только параллельных,
    но и последовательных, и их комбинаций:

    1 H2 + 2 NH3 + 1 N2 == 2 N2H4
    3 H2 + 1 N2 == 2 NH3
    1 H2 + 1 N2H4 == 2 NH3
    4 NH3 + 1 N2 == 3 N2H4
    2 H2 + 1 N2 == 1 N2H4

    ————
    2.
    В подразделе «4.1. Критерий равновесия — переход к дифференциалам» (с.27) задается уравнение (4.1) минимума энергии Гиббса равновесного состава гомогенной системы подобно следующему
    (к сожалению не могу здесь привести настоящую картинку формулы):

    min G = Sum_i ( n_i * mu_i )

    где mu_i = mu’_i + R*T*ln( n_i / n )

    Далее в подразделе «4.3. Пример: Расчёт равновесного состава при синтезе аммиака» (с.29) эта математическая задача (поиск условного минимума G) решается методом неопределённых множителей Лагранжа, а именно через частные производные по переменным (n_H2, n_N2, n_NH3, n_N2H4 и двум множителям Лагранжа — lambda_H, lambda_N ).

    Здесь получается система из шести нелинейных уравнений, например,
    для частной производной n_H2:

    df/dn_H2 = 0 = mu’_H2 + R*T*ln( n_H2 / n ) = 2 * lambda_H

    (аналогично и для других 4-х)

    Но для меня остается загадкой исчезновение слагаемых R*T,
    которые неминуемо должны были сопровождать первые 4-ре уравнения.

    Т.е. поясню, при нахождении частной производной, например, по n_H2, для:

    f = mu_H2 * n_H2 + mu_N2 * n_N2 + …
    + lambda_H * ( b_H — 2*n_H2 — … ) +
    + lambda_N * ( b_N — 2*n_N2 — … )

    Должно было получится:

    df/dn_H2 = mu’_H2 + R*T*ln( n_H2 / n ) — 2 * lambda_H + R*T

    поскольку mu_H2 зависит от n_H2 через логарифм.

    Поправьте меня пожалуйста. Может я в чём то ошибаюсь?

    ———

    Может я не совсем понял, но в работе http://evgenii.rudnyi.ru/doc/papers1/02pc_hg1201.pdf
    не обнаружил примера фазовой диаграммы Bi-Se.
    🙁

    Еще раз извиняюсь за беспокойство,
    Можно Вас попросить поделиться ссылкой на программу VCS
    (Villars-Cruise-Smith). Я буду Вам очень признательным.

    С assess я разбираюсь и пока не задаю лишних вопросов.

    Спасибо.
    С уважением,
    Александр Гороховский.

  2. Ссылка на Bi-Se была неправильная, я поправил. Спасибо.

    VCS из книги Smith&Missen:
    http://blog.rudnyi.ru/ru/2011/04/vcs-smith-missen.html
    У мнея к сожалению сейчас книги нет.

    Ссылка [5] уже морально устарела, нет никакого смысла ее искать.

    1) Вы правы можно составить множество реакций, не только две. Однако имеется в виду, что есть только две линейно независимые реакции, остальные — это их линейные комбинации. См. также
    http://blog.rudnyi.ru/ru/2011/04/khimicheskie-reaktsii-v-tdlib.htm

    2) Здесь еще сложнее, поскольку n_H2 такжe содержится в n:

    mu_i = mu’_i + R*T*ln( n_i / Sum_i n_i )

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

    P.S. Тема AmberPanther не показывает пустые строки в комментариях. Пока не знаю, как с этим бороться.

  3. Сейчас попробовал на http://www.wolframalpha.com/

    derivative (RT n1 log(n1/(n1 + n2)) + RT n2 log(n2/(n1 + n2)))

    Система выдает RT log(n1/(n1 + n2)). Так что должно быть все правильно.

  4. Да. Ваша правда. Недоглядел 🙂
    Проверил в своем любимом Maple всю систему:

    > f[H2] := n[H2] * ( mu[H2] + RT * ln( n[H2] / (n[H2] + n[N2] + n[NH3] + n[N2H4] ) ) );

    > f[N2] := n[N2] * ( mu[N2] + RT * ln( n[N2] / (n[H2] + n[N2] + n[NH3] + n[N2H4] ) ) ) ;

    > f[NH3] := n[NH3] * ( mu[NH3] + RT * ln( n[NH3] / (n[H2] + n[N2] + n[NH3] + n[N2H4] ) ) ) ;

    > f[N2H4] := n[N2H4] * ( mu[N2H4] + RT * ln( n[N2H4] / (n[H2] + n[N2] + n[NH3] + n[N2H4] ) ) ) ;

    > diff(f[H2] + f[N2] + f[NH3] + f[N2H4], n[H2]):
    > simplify(%);

    Получилась классика:
    mu[H2]+RT*ln(n[H2]/(n[H2]+n[N2]+n[NH3]+n[N2H4]))

    Отлегло.
    🙂

  5. У меня получилось рассчитать (с помощью assess) и построить (через gnuplot) фазовую диаграмму Bi-Se из каталога ex/bise. Нет слов — выглядит очень эффектно.

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

    Если не ошибаюсь, то в sys.mod задаются результаты аппроксимаций эксперимента в виде: состав — температура ?

    Или коэффициенты (A, B, C, …) этих полиномом вида
    A + B*T + C*T*ln(T) + D*T^2 + E*T^3 + F*T^7 + G/T + H/T^9

    получают из справочных данных?
    Может NIST ?
    (хотя там порядок полинома гораздо меньше)

  6. Я написал первое введение в TDLIB

    http://blog.rudnyi.ru/ru/2011/05/chistye-veshchestva.html

    Посмотрите пожалуйста. Там еще нет ответа на ваш вопрос, но это первый этап. Растворы будут в следующем введении.

  7. Вот бы еще понять, как в файле alg.mod задаётся описание возможных фазовых равновесий в системе. И на основании каких сведений сформировать представление о возможной будущей топологии диаграммы состояния, чтобы сформировать файл pd.out.mod ?

    А можно ли обойтись без задания топологии, чтобы TDLIB сам предложил возможные варианты?

  8. Я только что написал про растворы

    http://blog.rudnyi.ru/ru/2011/05/tdlib-regulyarnyi-rastvor.html

    Объект PhaseEquilibrium будет в следующем документе. К сожалению, топологию должен задавать пользователь, поскольку в TDLIB выполнено только решение уравнений 2.26.

  9. Премного благодарен.
    Возможно это будет выглядеть не в общей хронологии обсуждения,
    но чтобы не кануло в лету, а осталось для других кто пойдёт следом, приведу последовательность компиляции tdlib, которую Вы помогли осуществить:
    $ cd ./tdlib/lib/
    $ wget http://www.netlib.org/lapack/lapack-3.3.0.tgz
    $ tar zxvf lapack-3.3.0.tgz
    $ cd ./lapack-3.3.0

    $ cp INSTALL/make.inc.gfortran make.inc
    $ make blaslib
    $ make lapacklib
    $ ls *.a
    blas_LINUX.a lapack_LINUX.a

    $ cp blas_LINUX.a ../../bin/libblas.a
    $ cp lapack_LINUX.a ../../bin/liblapack.a
    $ cd ..

    Далее использовал gfortran, посему

    $ wget http://www.netlib.org/f2c/libf2c.zip
    $ mkdir libf2c
    $ unzip libf2c.zip -d ./libf2c
    $ cd libf2c
    $ cp makefile.u makefile
    $ make
    $ make hadd
    $ cp libf2c.a ../../bin/
    $ cd ..

    И наконец:

    $ make

    В итоге получил в ./bin/ долгожданный

    ./assess