Содержание: Термодинамика
http://evgenii.rudnyi.ru/doc/misc/tdlib.html
Файлы доступны в архиве http://evgenii.rudnyi.ru/soft/tdlib00+.tar.gz.
ПО TDLIB было написано мной в для решения разнообразных задач химической термодинамики, в том числе и расчете фазовых диаграмм. Однако ПО создавалось в ходе научной работе по так называемой оптимизации фазовых диаграмм, когда целью было определение энергий Гиббса из доступных экспериментальных данных, включая точки на фазовой диаграмме. Это наложило свой отпечаток на ПО и самое главное на приоритеты. В первую очередь программировались задачи, которые нужны были для решения текущей научной работы.
При компиляции библиотеки получается консольное приложение assess.exe, которое как раз может решать обратную задачу: определение энергии Гиббса из доступных экспериментальных данных. Однако, поскольку решение прямой задачи (расчет равновесного состава) является частью обратной задачи, то данное приложение можно использовать только для прямой задачи. Решение прямой задачи будет коротко рассмотрено в настоящем документе.
Должен отметить, что при расчете фазовой диаграммы наиболее оптимально использовать минимизацию энергии Гиббса, например см. ПО для расчета равновесного состава в трехкомпонентной системе 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), а также другие разделы).
Чистые вещества
Простейший объект в химической термодинамике — это стехиометрические фазы, которые не являются растворами. На жаргоне вычислительной термодинамики их также называют точечными фазами. Мольная энергия Гиббса в этом случае зависит только от температуры и давления G(T, p) и первой задачей TDLIB было дать возможность описывать произвольные функции от температуры и давления. Когда мольная энергия Гиббса известна, остальные термодинамические свойства можно рассчитать по уравнениям 2.9 (см. документацию (doc/00tdlib.pdf)).
Объект в TDLIB, моделирующий точечную фазу называется соответственно PointPhase. Все вычисления однако проводятся в объекте species и объект PointPhase является простым контейнером, который содержит один объект species. Причина для такого решения заключается в том, что при моделировании растворов требуются энергии Гиббса чистых компонентов, которые также выражаются как G(T, p). Введение специального объекта species позволило таким образом уменьшить количество необходимых объектов в библиотеке.
В простейшом виде PointPhase выглядит следующим образом
<PointPhase id=test>
<species> Y2O3
</species>
</PointPhase>
где объект species содержит химическую формулу вещества. Данный пример находится в директории ex/phase/pp/ в файле ex/phase/pp/ex1.mod. Можно посмотреть на свойства этого вещества посредством команды
$ assess ex1 out
T G H S Cp
300 0 0 0 0
600 0 0 0 0
900 0 0 0 0
1000 0 0 0 0
T V dVdT dVdp
300 0 0 0
600 0 0 0
900 0 0 0
1000 0 0 0
которая при помощи файла ex/phase/pp/out.mod распечатывает температурную таблицу свойств точечной фазы из первого файла. Как и следовало ожидать, все значению равняются нулю, поскольку мы еще не задали мольную энергию Гиббса, которая по умолчанию равняется нулю.
В общем случае для задания G(T, p) можно использовать объект calc_Tp, который является интерпретатором и позволяет задать произвольную функцию в переменных T и p. Пример в файле calc_tp.mod
<PointPhase id=test>
<species>
<calc_Tp>
(
10 + 5*T - 2*T*log(T) - 0.4*sqrt(T)) + R*T*log(p)
)
</calc_Tp>
</species>
</PointPhase>
Теперь
$ assess calc_tp out
T G H S Cp
300 -1919.2 606.536 8.41911 1.99423
600 -4676.11 1205.1 9.80202 1.99592
900 -7746.31 1804 10.6115 1.99667
1000 -8818.16 2003.68 10.8218 1.99684
T V dVdT dVdp
300 2494.32 8.31441 -2494.32
600 4988.65 8.31441 -4988.65
900 7482.97 8.31441 -7482.97
1000 8314.41 8.31441 -8314.41
Использование интерпретатора возможно, однако по разным причинам я счел это недостаточным (см. раздел 2.4 Some Thoughts on Possible Implementations в документации). В результате было создано много (наверное слишком много) объектов типа func_tp, перечисленных в таблице 3.2 на стр 26 в документации. Они более эффективны, чем calc_Tp, и их примеры можно найти в директории ex/phase/pp/. Примеры можно использовать аналогичным способом.
Для комбинирования объектов типа между собой можно использовать объекты complex_Tp и compound_Tp. Первый суммирует две объекта типа func_tp, второй позволяет составить составную функцию, которая в разных интервалах будет применять разные объекты. Здесь полезно выполнить следующие команды
assess cp_bb2 out
assess ideal_gas out
assess complex_tp out
assess compound_tp out
и взглянуть на соответствующие файлы.
Больше информации можно найти в разделах 3.2.3 — 3.2.5 на стр. 25-30 в документации. Раздел 3.2.5 можно пропустить, поскольку эта возможность была сделана для специальных приложений.
Необходимые численные данные для функций энергий Гиббса чистых веществ можно найти в термодинамических справочниках и термодинамических базах данных.
Регулярный раствор
В общем случае мольная энергия Гиббса является функцией концентрации раствора и в TDLIB есть несколько объектов для моделирования функции G(T, p, x1, …, xN), где в дополнение к температуре и давлению появляются мольные доли. Наиболее простая модель раствора описывает энергию Гиббса в виде уравнения 3.2 в документации (doc/00tdlib.pdf), где дополнительные члены, описывающую избыточную энергию Гиббса выражаются в виде полиномов, что приводит к названию полиномиальный раствор. В разделах 3.2.6 и 3.2.8 документации (стр 30-34) описывается объект SimpleSolution, который позволяет моделировать разные полиномы (Redlich-Kister, Borelius, Hoch-Arpshofen) и разные проекции (Muggianu, Kohler) в случае многокомпонентных растворов. Примеры, описанные в этих разделах, пожалуй несколько сложны, и ниже я опишу использование TDLIB в случае простейшего регулярного раствора (см. файлы в ex/regular_ab/L).
G = xA GA + xB GB + RT xA ln xA + RT xB ln xB + A0 xA xB
Файл ex/regular_ab/L/L.mod содержит соответствующий объект SimpleSolution, где параметр взаимодействия определен в начале файла. В файле энергии Гиббса чистых компонентов приняты за нуль и таким образом энергия Гиббса равна энергии Гиббса смешения. Вначале параметр взаимодействия равен нулю и мы имеем идеальный раствор. Команда
$ assess L out
показывает энергию Гиббса как функцию мольной доли при 1000 К.
Задачей будет построение энергии Гиббса при 1000 К, когда параметр взаимодействия меняется от -20 кДж до 20 кДж. Для этого придется вручную изменять величину параметра A в файле ex/regular_ab/L/L.mod в первой строке. Вначале изменим ее на -20000, сохраним изменения и подадим команду
$ assess L out -o g-20
которая сохранит энергию Гиббса в этом случае в файл g-20.G. Затем изменим величину параметра взаимодействия на -10000 и командой
$ assess L out -o g-10
сохраним энергию Гиббса в файл g-10.G. Повторим то же самое для 0 (файл g0.G), 10000 (файл g10.G) и 20000 (файл g20.G).
Теперь файл ex/regular_ab/L/plot содержит команду для gnuplot, которая покажет все энергии Гиббса вместе. Если все прошло удачно, то команда
$ gnuplot plot -
покажет вам результаты (см. ex/regular_ab/L/fig.png).
Можно увидеть, что большие положительные значения параметра взаимодействия меняют форму энергии Гиббса смешения — она становится невыпуклой. Это соответствует расслаиванию.
Энергии Гиббса растворов можно найти в статьях, где коэффициенты в модели раствора были найдены из доступных экспериментальных данных, или в термодинамических базах данных растворах.
Объект фазовые равновесия
Когда мольные энергии Гиббса фаз готовы, можно рассчитать фазовые равновесия посредством объекта PhaseEquilibrium. В документации (doc/00tdlib.pdf) ему посвящен раздел 3.4.7 на стр. 46-53. В настоящем разделе я кратко опишу объект PhaseEquilibrium на примере расчета ликвидуса в двухкомпонентном системе А-В.
В качестве примера я использую статью
Dmitri V. Malakhov, Thevika Balakumar
Undulate phase boundaries on binary T-x diagrams
Computer Coupling of Phase Diagrams and Thermochemistry 32 (2008) 89-93
В статье обсуждаются возможность существования точки перегиба на кривой ликвидуса и показывается, что точка перегиба возможна даже в случае идеального раствора. Ниже я покажу как с TDLIB можно построить рис. 2 из статьи.
С точки зрения математики объект PhaseEquilibrium позволяет описать систему нелинейных уравнений 2.26 (см. документацию) и решить ее. Файлы к примеру находятся в ex/regular_ab/malakhov08.
Мольные энергии Гиббса фаз описаны в ex/regular_ab/malakhov08/sys.mod, где L — это идеальный двухкомпонентный раствор A-B и A — это фаза А в кристаллическом состоянии. Предполагается, что твердые растворы отсутствуют и А — это точечная фаза. В файле задается два параметра, температура плавления и энтропия плавления. Переменная R=1 в глобальных переменных задает, что энтропия выражена в безразмерных единицах (S/R). Энтальпия плавления рассчитывается как Del_Hm = Tm Del_Sm.
Энергия Гиббса расплава А принята за нуль и таким образом энергия Гиббса твердого А задается как энергия Гиббса реакции затвердевания
А_l = A_s
которая в приближении Del_Cp = 0 выражается как
Del_G(T) = -(Del_Hm - T Del_Sm)
Система уравнений 2.26 (см. документацию) при постоянном давлении в нашем случае сводится к одному алгебраическому уравнению
Gm(As, T) = mu(Al, T, x2)
из которого при заданной мольной доли в расплаве можно вычислить температуру плавления (или из температуры плавления мольную долю расплава).
Объект PhaseEquilibrium находится в файле ex/regular_ab/malakhov08/alg.mod и выглядит следующим образом:
<PhaseEquilibrium id=A_L> <phases> <phase IDREF=L></phase> <phase IDREF=A></phase> </phases> <state status=constraint> x(L,B) </state> <state status=unknown> T </state> </PhaseEquilibrium>
Вначале идет список фаз, потом определяют зависимые и независимые величины. В данном случае, независимой величиной является мольная доля B, рассчитываемой величиной температура. В описании объекта ничего не сказано про давление и про мольную долю А, в отношении которых решения будут приняты по умолчанию. Команда
$ assess -m sys alg
позволяет вывести описание объекта в стандартном виде и узнать, какие решения по умолчанию были приняты. TDLIB проинтерпретировала ввод выше как
<PhaseEquilibrium id=A_L debug=0 SaveSolution=0 ThrowException=1> <phases> <phase IDREF=L></phase> <phase IDREF=A></phase> </phases> <state status=dependent> x(L,A) </state> <state status=constraint> x(L,B) </state> <state status=unknown> T </state> <state status=HardConstraint value=1> p </state> </PhaseEquilibrium>
Файл ex/regular_ab/malakhov08/out.mod организует цикл по мольной доле и выдает рассчитанные точки ликвидуса в формате gnuplot. Две команды
$ assess sys alg out -o 1
$ wgnuplot 1.pd -
позволяют построить рисунок ex/regular_ab/malakhov08/1.png, где показана линия ликвидуса в случае энтропии плавления, равной единице.
Следует отметить, что рассчитанная кривая включает в себя как стабильные так и метастабильные равновесия A-L. Например если мы включим в расчет твердый B, то ликвидус A-L ниже температуры эвтектики будет метастабильным по отношению к чистым твердым A и B.
Для построения семейства кривых требуется вручную менять значения энтропии плавления (см. выше раздел Регулярный раствор). С другой стороны, это также можно организовать путем показанным в ex/regular_ab/malakhov08/fig2.mod. Команды
assess fig2 -o out
gnuplot plot
дают рисунок в формате png (ex/regular_ab/malakhov08/fig2.png), который показывает, как линия ликвидуса зависит от энтропии плавления.
В документации в разделе 3.4.7 на стр. 46-53 можно найти другие примеры и дополнительные возможности объекта PhaseEquilibrium.
Уравнивание химических реакций
Уравнивание химических реакций, в особенности окислительно-восстановительных, это не самое приятное времяпровождение. В то же время существует способ, описанный в книге Smith & Missen (см. раздел ниже), сделать это совершенно формально на основе линейной алгебры. Коротко метод описан в
http://evgenii.rudnyi.ru/doc/ru/teaching/313/313group.pdf
на стр. 11 — 15. Формульная матрица (ex. 1 и ур. 2.1) содержит в себе всю необходимую информацию и стехиометрическую матрицу, которая содержит уравненные реакции, можно получить очень эффективно посредством LU-разложения. Алгоритм не обеспечивает получение целых стехиометрических коэффициентов. Если это необходимо, то надо можно сделать следующий шаг.
Естественно, что можно найти большее число реакций (любая линейная комбинация из столбцов стехиометрической матрицы также дает химическую реакцию), однако с практической точки зрения важен набор линейно независимых реакций. В термодинамике набор найденных реакций должен просто соответствовать материальному балансу и механизм как таковой не важен. С этой точки зрения наверное будет интересно посмотреть на набор ортонормированных реакций в
http://evgenii.rudnyi.ru/doc/papers1/93jct_va.pdf
(ур. 1-6) которые мне однажды потребовались из соображений математической статистики.
В TDLIB данный алгоритм использован в нескольких объектах. Наиболее наглядно это происходит в точечной фазе (PointPhase, см. стр. 29 в документации). Пример можно найти в ex/phase/pp/k2so4.mod. Команда
assess -m k2so4.mod
показывает уравненную реакцию. Если из списка компонентов уравнять реакцию не удается, то будет выдана ошибка. В объекте PointPhase ищется только одна реакция и лишние компоненты просто удаляются.
Далее алгоритм нахождения стехиометрической матрицы содержится в VCS, которая является частью объекта associated_solution, и в новом объекте AssociatedSolution. Объект PhaseEquilibrium также автоматически составляет уравнения между веществами в разных фазах. В данных объектах однако все происходит формально и скрыто от пользователя, некоторая информации содержится только в файлах отладки (debug=1).
Программа VCS из книги Smith&Missen
William R. Smith, Ronald W. Missen, Chemical Reaction Equilibrium Analysis: Theory and Algorithms.
Отличная книга, которая прекрасно описывает проведение расчета равновесного состава идеального ассоциированного раствора. На базе этой книги я сделал в свое время первую версию курса Математическое моделирование в химической термодинамике
http://evgenii.rudnyi.ru/doc/ru/teaching/313/313group.pdf
В книге также приведена программа VCS на языке Фортран, в которой запрограммирован алгоритм Villars-Cruise-Smith для проведения расчетов равновесного состава. Я очень долго использовал эту программу в своей работе (см. VCS в http://evgenii.rudnyi.ru/programming.html и соотвественно http://evgenii.rudnyi.ru/publications.html#plasma).
В TDLIB также была использована VCS (объект associated_solution), где вначале код был превращен в C посредством f2c, а затем был сделан класс C++. К сожалению, я не мог включить VCS в версию библиотеки, доступную в открытом доступе. Если вам потребуется объект associated_solution, напишите.