Физическая химия

Глава 9. Использование ПЭВМ при проведении термодинамических расчетов

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

9.1. Расчет параметров равновесного состояния гетерогенной системы на ЭВМ

Расчет параметров равновесного состояния природной гетерогенной системы рассмотрим на примере расчета равновесного состава грунтовых вод. Такой расчет может быть выполнен на основании моделирования процессов взаимодействия грунтовых вод с окружающими породами. В процессе установления равновесия происходит растворение одних веществ и выпадение в осадок других веществ. Одним из этапов анализа гетерогенного равновесия является решение задачи о формах нахождения элементов в природных водах для заданных параметров, например, при определенных значениях температуры и рН раствора и последующий анализ состава раствора при изменении этих параметров. В настоящий момент существуют два алгоритма решения данной задачи: 1) использование констант равновесия реакций, в результате которых устанавливается равновесие в системе; 2) расчет равновесного состава гидрогеохимической системы методом минимизации свободной энергии. Первый алгоритм реализован в программах под общим названием MIF (MIF1, MIF2 и т. д.). Второй алгоритм реализован в программах GIBBS, BALANS, MATRIX. Две последние программы разработаны на кафедре химии МГРИ (Н. Н. Акинфиев, В. А. Умрихин). В лекциях будет рассмотрен алгоритм расчета на основании минимизации потенциала Гиббса.

Постановка задачи. Пусть заданы внешние условия, в которых находится наша система, ее валовый состав и набор фаз, которые присутствуют в системе или которые могут возникнуть при эволюции системы. Требуется определить, какие конкретно фазы из заданного набора будут присутствовать в системе при равновесии. Например, пусть нас интересует состав насыщенного раствора кальцита в присутствии углекислого газа, который входит в состав атмосферы. При равновесии наша система будет трехфазной — насыщенный раствор кальцита в воде + твердая фаза (кальцит) + углекислый газ. Выполним расчет для стандартных условий: количество воды — 1 кг (55,51 моля), температура — 298 °К, давление в системе — 1 бар, причем парциальное давление углекислого газа 3 × 10–4 бара (среднестатистическое парциальное давление газа на поверхности Земли). Считаем, что раствор состоит из девяти компонентов:

Са2+, Са(ОН)+, СаСО03, Н2СО03, НСО3, СО32–, ОН, Н+, Н20.

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

Алгоритм решения задачи. Самопроизвольно протекающие химические реакции, приближающие нашу систему к состоянию равновесия, приводят к изменению компонентного состава системы и, следовательно, к изменению х термодинамических параметров системы. Та из термодинамических функций, которая достигает экстремального значения в состоянии равновесия, является термодинамическим потенциалом системы при заданных условиях, а задача нахождения экстремума термодинамической функции равносильна нахождению равновесного состава системы. Для изобарно-изотермических условий равновесный состав системы может быть найден путем минимизации энергии Гиббса. Таким образом, при заданных значениях давления, температуры и валового состава системы, система придет в состояние равновесия, когда энергия Гиббса системы достигнет своего минимального значения.

С математической точки зрения задача сводится к нахождению минимального значения функции многих переменных при заданных граничных условиях (валовом составе системы):

G = nj × µj min (9.1)

aij × mj = bi (9.2)

aij — количество независимого i-компонента в одном моле j-компонента системы;

mj — моляльная концентрация j-компонента;

bi — суммарное количество независимого i-компонента, равное его общей (валовой) концентрации.

Например, в задаче о растворении кальцита в присутствии углекислого газа, независимые компоненты (элементы, из которых можно создать любой компонент системы) — это кальций Са (i =1), углерод С (i = 2), водород Н (i = 3) и кислород О (i = 4); i — номер независимого компонента. j = 11 — число компонентов системы, aij — компонент матрицы баланса. Если j = 4, то это компонент раствора Н2СО30, его коэффициенты aij равны: а14 = 0, а24 = 1, а34 = 2, а44 = 3.

Задача нахождения минимума функции (9.1) при граничных условиях (9.2) эквивалентна нахождению минимума функционала (9.3):

(9.3)

Введем обозначение:

zj = µj/RT (9.4)

µj — химический потенциал;

j — компонента системы;

zj — приведенный потенциал компонента.

Таблица 9.1. Компонентный состав системы и матрица коэффициентов

Фазы (U)

Независимые компоненты (i)

(U = 1) раствор

1

2

3

4

j

(компоненты раствора)

Са

С

Н

О

1

Са+

1

0

0

0

2

Са(ОН)+

1

0

1

1

3

СаСО30

0

1

2

3

4

Н2СО30

0

1

2

3

5

НСО3

0

1

1

3

6

СО32–

0

1

0

3

7

ОН

0

0

1

1

8

Н+

0

0

1

0

9

Н2О

0

0

2

1

10 (U = 2) СаСО3 кальцит

1

1

0

3

11 (U = 3) СО2 углек. газ

0

1

0

2

Тогда из уравнения (9.5) можно вычислить концентрации ком­понентов раствора при равновесии, а из уравнения (9.6) приведенный химический потенциал компонента раствора.

mj = exp( – zj + Σui × aij) (9.5)

zj = Σui × aij (9.6)

В таблице 9.2 приведены химические потенциалы независимых компонентов ui (их еще называют начальные приближения) задачи о растворении кальцита, а в таблице 9.1 компонентный состав раствора.

Таблица 9.2. Начальные приближения (ui) для независимых компонентов.

Са

С

О

Н

–230,708

–52,249

–57,473

–19,102

В таблице 9.3 приведены результаты расчета поставленной задачи.

Таблица 9.3. Результаты расчета компонентного состава раствора (моль/кг)

Са2+

Са(ОН)+

СаСО30

Н2СО30

НСО3

СО32–

4,49 × 10–4

1,24 × 10–8

6,11 × 10–6

1,02 × 10–5

8,79 × 10–4

8,16 × 10–6

кальцита растворилось 5,44 × 10–4 (моль/кг)

углекислого газа растворилось 5,41 × 10–4 (моль/кг)

рН раствора 8,29

9.2. Расчет компонентного состава природных вод на основании элементного анализа пробы грунтовых вод

Большое значение в оценке пригодности природных вод для хозяйственной деятельности имеет компонентный состав конкретного водного бассейна. При исследованиях в области охраны окружающей среды первостепенное значение имеет то обстоятельство, что разные миграционные формы одного и того же элемента имеют резко различную токсичность. Например, неорганические соединения трехвалентного мышьяка As(III) на порядок более токсичны, чем соединения пятивалентного мышьяка As(V).

Формами миграции элемента будем называть сложные и простые частицы, в состав которых входит данный элемент и которые не могут существовать вне раствора. Например, формы миграции кальция в растворе — это ионы Са2+, неорганические комплексные соединенияCa(OH)+, СаНСО3+, СаСО30, CaSO40, Ca(SO4)22–, органические комплексные соединения, например, с фульвокислотой СаФК0 и т. д. Движение (миграция) в пространстве таких частиц происходит вместе с раствором. При изменении параметров раствора возможно выпадение в осадок твердой фазы, в состав которой входит данный элемент. Тогда говорят о геохимических барьерах миграции данного элемента. Часто таким барьером является температура и рН раствора.

Постановка задачи. Маломинерализованная вода четвертичного горизонта белорусской части Полесья имеет состав, мг/кг:

Na + K

Са

Mg

Fe

Cu

Mn

H2CO30

SO4

Cl

pH

13,0

10,0

6,06

1,1

0,01

1,3

78,0

3,0

13,0

6,5

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

Алгоритм решения задачи. Алгоритм решения задачи мы сформулировали при анализе гетерогенного равновесия (см. 9.2).В соответствии с этим алгоритмом записываем матрицу баланса (табл. 9.4). При записи матрицы баланса мы считаем, что в растворе натрий, калий и хлор находятся в виде простых ионов Na+, K+, Cl в количествах, которые были определены для соответствующих элементов. Так как раствор практически нейтральный (рН = 6,5) не имеет смысл учитывать частицы типа Са(ОН)20, концентрации которых в нейтральных растворах малы. В таблице 9.4 приведены результаты расчета компонентного состава раствора. (Расчет выполнен с использованием программы «GROUND», кафедра химии МГРИ, Умрихин В. А.) Чтобы оценить, находится ли наш раствор в равновесии относительно, например, доломита CaMg(CO3)2 надо записать реакцию диссоциации минерала и вычислить константу

Таблица 9.4. Матрица баланса и результат расчета компонентного состава раствора

Компонент

Ca

Cu

Mn

Mg

Pb

S

ФК

Fe

C

O

H

мг/кг

1. Ca2+

1

0

0

0

0

0

0

0

0

0

0

9,8

2. Cu2+

0

1

0

0

0

0

0

0

0

0

0

0,0046

3. Mg2+

0

0

0

1

0

0

0

0

0

0

0

5,9

4. Mn2+

0

0

1

0

0

0

0

0

0

0

0

1,3

5. Pb2+

0

0

0

0

1

0

0

0

0

0

0

0,00366

6. SO42–

0

0

0

0

0

1

0

0

0

4

0

3,01

7. ФК2–

0

0

0

0

0

0

1

0

0

0

0

2,01

8. Fe2+

0

0

0

0

0

0

0

1

0

0

0

1,1

9. Ca(SO4)0

1

0

0

0

0

1

0

0

0

4

0

0,15

10. Ca(HCO3)+

1

0

0

0

0

0

0

0

1

3

1

0,47

11. CaCO30

1

0

0

0

0

0

0

0

1

3

0

0,00646

12. CaФК0

1

0

0

0

0

0

1

0

0

0

0

3,1 × 10-24

13. Mg(CO3)0

0

0

0

1

0

0

0

0

1

3

0

0,00313

14. Mg(HCO3)+

0

0

0

1

0

0

0

0

1

3

1

0,27

15. MgФК0

0

1

0

0

0

0

1

0

0

0

0

1,4 × 10-24

16. FeФК0

0

0

0

0

0

0

1

1

0

0

0

3,3 × 10-25

17. Fe(OH)+

0

0

0

0

0

0

0

1

0

1

1

0,00129

18. Mn(SO4)0

0

0

1

0

0

1

0

0

0

4

0

0,0234

19. Mn(HCO3)+

0

0

1

0

0

0

0

0

1

3

1

0,0222

20. HSO4

0

0

0

0

0

0

1

0

0

4

1

6,27 × 10–5

21. Pb(OH)+

0

0

0

0

0

1

0

0

0

1

1

0,00686

22. CuCO30

0

1

0

0

0

0

0

0

1

3

0

0,0219

23. CO32–

0

0

0

0

0

0

0

0

1

3

0

0,0137

24. H2CO30

0

0

0

0

0

0

0

0

1

3

2

58

25. HCO3

0

0

0

0

0

0

0

0

1

3

1

78

26. OH

0

0

0

0

0

0

0

0

0

1

1

5,72 × 10–4

27. H+

0

0

0

0

0

0

0

0

0

0

1

3,60 × 10–4

28. H2O

0

0

0

0

0

0

0

0

0

1

2

1000

равновесия реакции: CaMg(CO3)2 = Ca2+ Mg2+ 2CO32–. Константа равновесия реакции K = [Ca2+] × [Mg2+] × [CO32]2, уравнение Вант-Гоффа:

Результат расчета: –lgК = 15,44. Отметим, что десятичный логарифм, взятый с обратным знаком, константы реакции диссоциации плохо растворимого вещества называется произведением растворимости ПР этого вещества. Таким образом, ПР(доломита) = 15,44. По результатам расчета компонентного состава раствора, вычислим ионное произведение для доломита: ИП(доломита) = [Ca2+] × [Mg2+] × [CO32–]2 = 2,45 × 10–4 × 2,46 × × 10–4 × (2,28 × 10–7)2 = 3,13 × 10–21. lg[ИП] = –20,5. Так как ИП много меньше, чем константа равновесия реакции диссоциации доломита (К = 3,63 × 10–15), мы можем утверждать, что в водах Полесья карбонаты будут растворяться. Чтобы оценить степень насыщенности раствора относительно атмосферной углекислоты, вычислим концентрацию частицы Н2СО30 при равновесии раствора с углекислым газом и сравним с фактической концентрацией этой частицы в растворе.

СО2газ + Н2О = Н2СО30

[H2CO30] = K × P(CO2)

Уравнение Вант-Гоффа:

парциальное давление углекислого газа в атмосфере Р(СО2) = 3 × 10–4 бара = 30 Па. Резуль­таты расчета: К = 0,0293; [H2CO30] = 0,0293 × 30 = 0,879 (моль/кг) = 54 мг/кг (концентрация углекислоты при равновесии). Воды Полесья практически находятся в равновесии с углекислотой [H2CO30]факт = 58 мг/кг (см. табл. 9.4).

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