История развития методов физико-химического моделирования

Глава 2. Физико-химическое моделирование и современные вычислительные средства

История развития методов физико-химического моделирования

В 60-70-х годах XX века благодаря исследованиям P.M. Гаррелса (Гаррелс, 1960, 1962) и Ч.Л. Крайста (Гаррелс, Крайст, 1968) в геохимии широкое распространение получил метод анализа полей устойчивости минералов и компонентов водного раствора в координатах независимых или условно-независимых параметров состояния системы (Говоров, 1970; Дроздовская, Мельник, 1967, 1968, 1975, 1986; Иванова, Ходаковский, 1968; Карпов, 1965; Колонин, Птицын, 1974; Коренбаум, 1970; Мельник, 1968; Browne, Ellis, 1970; Tardy, 1971). Построение таких диаграмм опирается на систему уравнений действия масс изучаемой модели. Система же уравнений баланса масс в целом здесь не учитывается. Последнее обстоятельство является главным недостатком указанного метода. Использование только системы уравнений действия масс не позволяет в большинстве случаев решать однозначно задачу количественного определения равновесного фазового и компонентного состава моделируемой системы по ее исходному химическому составу.

В исследованиях Р.П. Рафальского (Рафальский, 1973) и И.К. Карпова (Карпов, 1981) математически строго на примере равновесного состава систем S-Н-О и Н-N-О-С показано, что способ изучения физико-химических моделей в геохимии, ориентированный на одностороннее использование уравнений действия масс, не может стать основой физико-химического моделирования на ЭВМ. Ограниченные возможности «диаграммного метода» объективно вытекают из термодинамического содержания теоремы Дюгема (Пригожин, Дефей, 1966), смысл которой заключается в том, что полное количественное описание закрытых систем, в том числе их компонентного и фазового состава, невозможно получить посредством задания только интенсивных параметров. Сложившаяся практика построения диаграмм Eh-рН и других невольно приводит к перестановке причины и следствия.

Координатная система рассматривается уже как независимый фактор состояния минерального парагенезиса и состава, равновесного с минералами водного раствора и газовой фазы. На самом деле, в природных закрытых системах, условия равновесия в которых подчинены локальному принципу Д.С. Коржинского, не рН определяет минеральный состав парагенезиса, а, наоборот, химический состав горной породы при данных Р и Т или Т и V определяет ее минеральный парагенезис, газовую фазу и состав равновесного с ним водного раствора, значения рН, а также активности других компонентов, фигурирующих на диаграммах в качестве координатных переменных.

Самым простым методом расчета химического равновесия, в котором одновременно учитывается уравнение действия и баланса масс, является метод сокращения числа неизвестных путем последовательных подстановок. Таким образом, получают одно – два нелинейных уравнения, иногда весьма сложных и громоздких (Хельгесон, 1967), с соответствующим числом неизвестных. Эти итоговые уравнения решают различными численными и графическими методами, а затем обратными подстановками находят остальные неизвестные. Большое значение в этом методе имеют физические соображения и эмпирический опыт исследователя. Например, при решении нелинейного уравнения с одним неизвестным часто возникает проблема выбора из нескольких решений, математически удовлетворяющих этому уравнению, одного, соответствующего ожидаемому результату. Примеры применения этого метода изложены в специальных руководствах и учебниках (Батлер, 1973; Введенский, 1960; Карапетьянц, 1953; Путилов, 1971; Garrels, Thompson, 1962, Suess, 1962).

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

В геохимии и петрологии метод подстановок успешно использовался разными исследователями. Г. Хельгесон (Хельгесон, 1967) исследовал комплексообразование в гидротермальном флюиде. P.M. Гаррелс и Ч.Л. Крайст (Гаррелс, Крайст, 1968) рассчитали компонентный состав модели морской воды и соленого озера. P.M. Гаррелс (Гаррелс, 1962) широко использовал метод подстановок для изучения разнообразных моделей равновесия карбонатов. Р.П. Рафальский (Рафальский, 1973) произвел тщательную оценку состава некоторых гидротермальных систем.

С.Р. Бринкли (Brinkley, 1946, 1947) была разработана универсальная аналитическая процедура, которая сводит постановку и численное решение задач химической термодинамики с учетом уравнений действия и баланса масс к выполнению определенной последовательности стандартных операций. Таким образом, удается избежать отдельных недостатков метода подстановок. Это достигается тем, что все компоненты системы подразделены на зависимые и независимые компоненты. Важно то, что сама процедура разделения компонентов у С.Р. Бринкли также формализована. Благодаря этому ему полностью удалось решить сложную задачу аналитического построения исходной системы стехиометрических уравнений реакций независимо от конкретных особенностей физико-химической модели. В отличие от метода подстановок метод Бринкли дает общий аналитический прием построения и численного решения различных задач физико-химического моделирования, т. е. он является универсальным. С.Р. Бринкли впервые ввел и практически продемонстрировал принцип алгоритмического решения задач химической термодинамики, что является одним из основных условий эффективного применения ЭВМ в изучении сложных химических равновесий. В дальнейшем метод Бринкли был усовершенствован и модифицирован (Бугаевский, Мухина, 1976; Круглов, Бугаевский, 1975, 1976).

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

В 1958 г. В. Вайтом с соавторами (White, Johnson, Dantzig, 1958) была впервые показана практическая возможность численной минимизации свободной энергии для расчета химического равновесия идеальной газовой смеси. В отличие от расчетной схемы С.Р. Бринкли и ее модификаций в методе Вайта становится ненужной предварительная работа по составлению стехиометрических уравнений реакций, выбору среди них линейно-независимой базы и расчету констант равновесия. Такой подход в значительной степени упрощает расчет сложных химических равновесий (Лузанов, 1966; Kandiner, Brinkley, 1950). Поэтому понятен тот резонанс, который вызвала статья В. Вайта (White, Johnson, Dantzig, 1958) среди физико-химиков, работающих в области теории и практики расчетов сложных химических равновесий (Dantzig, Dehaven, 1962; Naphtali, 1959, 1961; Raju, Krishnaswami, 1956). Вскоре метод минимизации был распространен на изучение химических равновесий при разнообразном выборе пар независимых параметров состояния системы (Zeleznik, Gordon, 1968).

В науках о Земле метод минимизации был впервые применен Е.Ф. Хильдом и Дж. Найфтоном (Heald, Naughton, 1962) для расчета химического равновесия в гомогенных смесях вулканических газов. И. Шимазу (Shimazu, 1967) использовал метод минимизации для расчета различных вариантов первичного состава планетных образований Солнечной системы. Во всех этих приложениях использовалась вычислительная процедура наискорейшего спуска, предложенная В. Вайтом (White, Johnson, Dantzig, 1958) первоначально для расчета равновесного состава гомогенной газовой фазы и обобщенная затем на гетерогенные системы Ф. Бойтоном (Boynton, 1960).

Более полное представление о существующих методах численного решения задач химического равновесия, в том числе метода минимизации, можно получить, обратившись к специальным работам (Волков, Рузайкин, 1974; Алемасов и др., 1971).

Итак, казалось бы, развитие методов термодинамического изучения физико-химических систем открыло принципиально новые возможности. Стратегическое направление, в котором ученые должны совершенствоваться, определилось. Первоначально – от метода диаграмм (Гаррелс и Крайст, 1968) – схем с фиксированным химическим составом PТ-условиями к аппарату стехиометрических реакций, опирающемуся на систему уравнений действия масс (Хельгесон, 1967; Garrels, Thompson, 1962). Затем – к универсальной аналитической процедуре построения и численного решения задач физико-химического моделирования с учетом уравнений действия и баланса масс (Brinkley, 1946, 1947). А в конечном итоге – отказ от громоздкого аппарата стехиометрических уравнений реакций, и переход к методу минимизации свободной энергии (White, Johnson, Dantzig, 1958), позволяющему представлять физико-химические модели в терминах теории математического программирования с учетом всех требований химической термодинамики (Карпов, 1981).

Тем не менее, особенности природных геохимических процессов, налагающие специфические ограничения на химическую термодинамику, по крайней мере, до основополагающих работ Д.С. Коржинского, не рассматривались. И только в работах И.К. Карпова было показано, что метод расчета химических равновесий путем численной минимизации свободной энергии системы может быть, с учетом всех требований химической термодинамики природных процессов, обобщен в стройную теорию приложения математического программирования к широкому классу задач (Карпов, 1965, 1968, 1971, 1977, 1981; Карпов, Киселёв, Дорогокупец, 1976; Карпов, Киселев, Летников, 1971, 1976; Карпов, Трошина, 1967, Карпов и др. 1973, 1974, 1991, 1998).

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

В настоящее время оба эти направления продолжают равноправно развиваться как у нас в стране, так и за рубежом. Воспользуемся обзорным описанием компьютерных программ (табл. 1), используемых в задачах взаимодействия «вода – горные породы», представленных С.Р. Крайновым (Крайнов, 1993). Обзор объективно свидетельствует о том, что компьютерные программы расчета по реакциям узко специализированы как по возможности изучения различных процессов, так и по функциональным характеристикам. Дело в том, что каждая программа расчета равновесия по реакциям имеет фиксированный реестр веществ, участвующих в расчете равновесия. Такой подход к программному обеспечению приводит к жесткой зависимости результативности моделирования от возможностей конкретной специализированной программы. Уже не слишком принципиальные изменения в постановке вызывают необходимость модификации программы или обращение к другой специализированной программе.

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

Принцип «каждой задаче своя программа» оправдан, когда ставится и решается сравнительно простая однотипная задача. Но он неприемлем, если нужно имитировать сложные физико-химические процессы, протекающие в реальных геологических условиях или в геотехнических системах.

Возможности метода минимизации успешно использовались в исследованиях, проводившихся на кафедре геохимии МГУ, в решении проблемы термодинамического моделирования и построения на этой основе модели (Гричук, 2000). Термодинамические расчеты выполнялись с помощью программного комплекса (ПК) GBFLOW, ориентированного на расчет равновесно-динамических моделей. ПК включает блок минимизации свободной энергии системы, разработанный Ю.В. Шваровым (Шваров, 1999), и блок изотопно-химического моделирования.

Таблица 2.1