АЛГОРИТМ РЕШЕНИЯ ОБРАТНЫХ ЗАДАЧ ПРИ ИЗУЧЕНИИ СКОРОСТИ РАСТВОРЕНИЯ МИНЕРАЛОВ В ВОДНЫХ РАСТВОРЫХ

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

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

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

Методика решения подобных задач (прямой и обратной) была опробована при изучении природных процессов: при обработке данных откачек из однородного пласта с покровом и в двухслойном потоке около реки

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

К построению алгоритма решения можно подойти со статистических позиций. Пусть y1,y2,y3...yn – некоторая заданная последовательность чисел, например, экспериментальные данные по временной или пространственной координате. Если наблюдения yi проводятся в точках пространства и (или) времени xi = (xi1,...,xik), то каждое наблюдение можно представить как сумму регулярной детерминированной составляющей φ( (xi1,...,xik,α) и случайной помехи:

(1)

где α – неизвестный многомерный параметр,εi – случайные ошибки, о которых можно предположить, что они имеют нулевое математическое ожидание, постоянную дисперсию и не коррелируют друг с другом. Тогда φ можно рассматривать как функцию регрессии. Если она линейна относительно неизвестного параметра α, то уравнение (1) называют моделью линейной регрессии; если же φ нелинейно относительно неизвестного параметра α, то – (1) является моделью нелинейной регрессии. Задача состоит в оценивании неизвестного параметра α. В качестве метода оценивания, как правило, применяется метод наименьших квадратов. При этом значение α, при котором достигается минимум функции качества F(α) = ∑(yi – φ(xi,α))2, называется оценкой метода наименьших квадратов (МНК). Если φ – линейная функция, то полученная оценка параметра α является несмещенной и имеющей наименьшую дисперсию. Если же известно, что случайная ошибка ε распределена нормально, то оценка МНК совпадает с оценкой неизвестного параметра α по методу наибольшего правдоподобия и для нелинейной функции φ.

Нелинейная задача о наименьших квадратах требует особых методов решения, Этот класс задач относится к области математики – обратные задачи. Под прямой задачей понимается аналитическое, численно-аналитическое, конечно-разностное (или в конечных элементах) решение задачи, например, гидродинамики подземных вод, сорбции, т.е. для заданного параметра α в точках наблюдений вычисляется функция φ(xi,α). Под обратной задачей понимаем поиск такого параметра α, при котором достигается минимум функции качества F(α).

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

Примером метода минимизации является метод Гаусса – Ньютона (или его модификация – метод Левенберга–Марквардта)

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

Использована термостатированная ячейка (V = 89 мл, 25°C), на дне которой размещали образец породы (S = 9.62 см2). Раствор перемешивали мешалкой. Для регистрации кинетических кривых применяли метод кондуктометрии. Образцы для опытов: P1ir, плотная порода бледно-голубого цвета с гетерогенно-блоковой текстурой, содержание гипса и ангидрита G = 7-78 и A = до 95 масс.% соответственно. Более подробно установка и методика проведения опытов описаны в работе

Скорость процесса растворения опытных образцов (R) рассчитывалась по выходу в раствор ионов Са2+ (С), как с одной, так и двумя переменными соответственно, модель 1 и – 2. Первая модель составлена на основе балансового уравнения кинетики одновременного растворения гипса и ангидрита (R = k1(Cm1 – C) + k2(Cm2 – C)2 , где k – константа скорости растворения, Cm – концентрация насыщения, индексы 1–2 характеризуют величины Cm, k, при растворении гипса и ангидрита, соответственно). При определенных заменах переменных и подстановках, получено общее уравнение Риккати с известными аналитическими решениями

Вторая модель – система двух дифференциальных уравнений, описывающих раздельное растворение гипса (R1 = k1(Cm1 – C1)) и ангидрита (R2 = k2(Cm2 – C2)2). В этих уравнениях переменные связаны функциональными зависимостями (например, Cm1(C2)), а содержание ионов Са2+ в растворе (С, С = С1 + С2) характеризует собой суммарное количество массы растворенных минералов

Обработка результатов опытов. При решении прямой задачи величины С, С1 и С2 определялись при заданных значениях параметров – k1, k2, Cm1, Cm2. При решении обратной задачи – в качестве результата принято значение k2, соответствующее минимуму функции качества F(α). Величина F(α) определялась, как сумма квадратов разности между модельными и измеренными концентрациями по всем расчетным моментам времени. Значения величины F(α) находились путем многократных решений прямой задачи при изменении начальных величин искомых параметров.

Для расчетов С, С1, С2 и k2 использовались авторские алгоритм и программа

Рисунок 1 - Зависимость величины R от C и минерального состава образцов (2–5)

Примечание: линии и маркера - значения, рассчитанные по моделям 1 и 2, соответственно: 2 = G – 7 и A – 90, 3 = 31 и 69, 4 = 74 и 26, 5 = 78 и 22 масс.%

DOI:10.60797/CHEM.2024.2.1.1

– Обработка опытных данных в виде зависимостей С(t) показала, что использование моделей 1 и 2 для расчета кинетики растворения гипсоангидритов в воде вполне оправданно (25 °C; С <0.01 ммоль/см3, (рис. 1).

– Рассчитанные значения величины C наиболее близки экспериментальным данным с учетом значений: k1 == 9.64 · 10–4 см/с; Cm1 = 0.0151 – 0.1598 · С2, ммоль/см3; Cm2 = 0.0194 ммоль/см3; r1 = 1; r2 = 2. При 25 °C k2 = 0.00095 см4/(ммоль · с) (рис. 1).

– Средние значения величины k2, рассчитанные по этим двум моделям различаются не более, чем на 40 %. Точность определения параметров подобных k2 ≈ ± (10-30%), поэтому можно считать, что обе предложенные модели применимы для интерпретации процессов совместного растворения гипса и ангидрита.

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

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

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

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