CC BY
9Математические заметки СВФУ Апрель—июнь, 2023. Том 30, № 2
УДК 519.63
ИДЕНТИФИКАЦИЯ СКОРОСТЕЙ ГОМОГЕННО-ГЕТЕРОГЕННОЙ РЕАКЦИИ В МАСШТАБЕ ПОР В ПОРИСТЫХ СРЕДАХ В. В. Григорьев
Аннотация. Представлена модель гомогенно-гетерогенной реакции в масштабе пор, основанная на уравнениях Стокса и уравнениях конвекции-диффузии-реакции с граничными условиями третьего рода на границах включений. Гомогенная реакция описывается как кубический автокатализ на всем поровом пространстве, а кинетика гетерогенной реакции описывается изотермой Ленгмюра. Численное решение задачи производится методом конечных элементов на кусочно-линейных элементах. Для дискретизации по времени используется схема Кранка — Николсон. Нелинейная задача решается итерационным методом Ньютона. Массоперенос смоделирован с рассчитанным полем скорости. Дополнительно проведен анализ чувствительности модели к параметрам для изучения их влияния на реагирующий перенос через пористую среду. Представлено численное решение обратной задачи, а именно, идентификация ключевых параметров, характеризующих реагирующий перенос на основе двух кривых проскока двух разных растворов. Рассмотрены зашумленные данные с разными амплитудами шума, включая смешанные амплитуды. Для приближенного решения многомерной обратной задачи применен метаэвристический Алгоритм Искусственной Пчелиной Колонии, который показал хорошую эффективность при достаточно малых вычислительных затратах.
Б01: 10.25587/8УРи.2023.74.45.008 Ключевые слова: гомогенно-гетерогенная реакция, пористые среды, масштаб пор, идентификация параметров, метод конечных элементов.
1. Введение
Реагирующий перенос в пористых средах встречается во многих промышленных задачах, например: сбор и захоронение С02, перенос загрязняющих веществ, каталитические фильтры, кислотное воздействие на нефтяные пласты, очистка сточных вод и т. д. Во всех этих перечисленных задачах присутствуют сложные химические процессы, такие как перенос реагентов и продуктов реакции, взаимодействие между разными жидкостями и химические реакции на границах раздела фаз [1,2]. Различают два вида реакций в реагирующих потоках: гетерогенный и гомогенный. Гомогенная реакция - реагирующие компоненты находятся в одной фазе, реакция происходит по всему объему. Гетерогенная реакция - реагирующие компоненты находятся в разных фазах, реакция
Работа выполнена при поддержке гранта Главы Республики Саха (Якутия) (Соглашение № 10 от 17 мая 2022 года) молодым ученым и гранта Российского Правительства, нацеленного на поддержку молодежных лабораторий (Р8И.0-2021-0015).
© 2023 Григорьев В. В.
происходит на границе раздела фаз [3]. Реагирующие процессы, как правило, происходят в масштабе пор, поэтому необходимо проводить исследования в этом масштабе для качественного и количественного понимания сложных процессов реакции.
Существует несколько способов провести экспериментальные исследования кинетики реакций в объекте исследования, некоторые из них описаны в [4,5]. В общей схеме через объект исследования пропускают раствор или газ с известной концентрацией, попутно фиксируя концентрацию на выходе. Эта динамика в последующем может сыграть роль дополнительной информации для решения задачи идентификации ключевых параметров кинетики реакции. Для получения надежных кинетических данных реакций необходимо очень точно контролировать условия реакции в экспериментальных установках. Это сложно по трем основным причинам: во-первых, одновременно с рассматриваемой химической реакцией происходят различные другие процессы, которые мешают контролю условий реакции; во-вторых, концентрация и давление не могут быть измерены непосредственно в образце; в-третьих, условия реакции могут происходить неравномерно внутри образца из-за морфологических неоднород-ностей [5]. Разработка эффективных вычислительных алгоритмов идентификации ключевых параметров, определяющих кинетику реакций, которые учитывали бы вышеизложенные нюансы экспериментальных исследований, имеет несомненную актуальность.
Численное моделирование как метод описания и изучения реагирующего переноса в пористых средах широко применяется в последние два десятилетия [6]. В целом численное моделирование процессов переноса в пористых средах можно классифицировать двумя классами: масштаб континуума и масштаб пор. В моделях масштаба континуума пористая среда обычно рассматривается как однородная и изотропная область, при этом неоднородностью на уровне пор пренебрегают. Решаются макроскопические уравнения сохранения, которые получаются на основе теорий объемного усреднения при задании феноменологических параметров для неявного учета микроструктуры пористой среды. Для однофазного потока жидкости были приняты уравнение Дарси, Бринкмана с расширением Дарси, Форхгеймера с расширением Дарси. Хотя эффективные транспортные свойства и их крупномасштабная пространственная неоднородность часто фигурируют в качестве важных факторов, определяющих транспортные процессы, признано, что микромасштабная неоднородность также играет важную роль, а процессы на уровне пор сильно влияют на производительность системы, долговечность и стоимость [7, 8].
При моделировании в масштабе пор микроскопические пористые структуры пористой среды определяются в явном виде, что обеспечивает детализацию распределения важных переменных на уровне пор, позволяет напрямую связать сложные процессы переноса с реалистичными пористыми структурами и, таким образом, могут обеспечить глубокое понимание взаимосвязей между структурами, процессами и характеристиками. Детальное моделирование в масштабе пор
требует эффективного метода трехмерной реконструкции. Для реконструкции одним из основных препятствий является низкое разрешение трехмерной микроструктуры. Несмотря на все большую доступность устройств для компьютерной томографии, их разрешения все еще недостаточно [9]. Прямую визуализацию трехмерной микроструктуры пористой среды получают с помощью экспериментальных методов визуализации, таких как электронная томография, ХСТ и ИБ-БЕМ. Пористая среда визуализируется с разных направлений для создания последовательных двухмерных поперечных сечений, которые впоследствии объединяются вместе для формирования трехмерной структуры. Описанные методы помогают хорошо понять морфологию пористых структур, однако для оценки изменчивости, связанной с геометрией и составляющими пористых сред, требуются многочисленные эксперименты для анализа статистических характеристик, которые все еще остаются весьма дорогим и времязатратным процессом [10].
Данная работа связана с решением прямой задачи и решением обратной коэффициентной задачи. Преследуются две цели: реализация математической модели, описывающей процесс гомогенно-гетерогенной реакции, и вычислительная идентификация скоростей гомогенно-гетерогенной реакции в масштабе пор. Процесс рассмотрен как изотермический, гетерогенная реакция описывается изотермой Ленгмюра [11], гомогенная реакция представляется как кубический автокатализ [12]. Ранее другие авторы уже решали гомогенно-гетерогенную реакцию в пористых средах: в работе [13] авторы исследовали течение наножид-кости в пористой среде; в [14] авторы рассмотрели магнитогидродинамический поток жидкости Кассона; в [15] исследован теплообмен в ферромагнитном потоке с химическими реакциями. Во всех перечисленных работах присутствует кубический автокатализ и реакция первого порядка, описывающая гетерогенную реакцию, которая состоит только из одного управляющего параметра. В нашей работе гетерогенная реакция описывается изотермой Ленгмюра, которая является более сложной трехпараметрической моделью и применима для описания широкого спектра адсорбентов [16]. Ранее нами были выполнены несколько работ по идентификации параметров скоростей реакций [17,18], также есть работы и по изучению отдельных факторов, влияющих на гетерогенную реакцию [19]. Во всех этих работах идентификация проводилась только для гетерогенных реакций. В этой работе идентификация будет проводиться сразу для двух типов реакций, происходящих одновременно.
Численное решение проводится при конечно-элементной аппроксимации по пространству. Поток жидкости описывается стационарными уравнениями Стокса, а реагирующий перенос описывается уравнением конвекции-диффузии-реакции. Для гидродинамических процессов используются элементы Тейлора — Худа (Р2 — Р1), а для переноса реагирующих веществ — полиномы Лагранжа первой степени (Р1). Дискретизация по времени проводится симметричной схемой Кранка — Николсон, которая абсолютно устойчива для линейных задач и имеет второй порядок точности. Нелинейная задача решается итерационным
методом Ньютона. Программный комплекс написан с использованием вычислительной платформы FEniCS [20]. Визуализация полученных результатов проводится с помощью программы Paraview и библиотеки matplotlib. Вычисления представлены для двумерной задачи (сечение), саму математическую модель и вычислительный алгоритм практически без изменений можно применить и для полноценного трехмерного случая.
Решение обратной коэффициентной задачи предложено искать с помощью метаэвристического метода. Метаэвристика — это не зависящая от задачи алгоритмическая структура высокого уровня, предназначенная для поиска решения оптимизационной задачи при ограниченных вычислительных мощностях с неполной или несовершенной начальной информацией. Существуют довольно много метаэвристических алгоритмов и всевозможные их модификации. Дополнительно с этим классом алгоритмов можно ознакомиться в [21, 22]. Будет использоваться Алгоритм Искусственной Пчелиной Колонии (АИПК, англ. Artificial Bee Colony Algorithm) [23], который хорошо себя показал в промышленных и инженерных задачах. Ранее в работе по применению метаэвристических алгоритмов была предложена модификация Пчелиного Алгоритма [24], который оказался эффективен, когда необходимо идентифицировать только два или три параметра. Для четырех и более параметров он неэффективен, поскольку делает много вызовов минимизируемой функции из-за отсутствия глобального поиска. АИПК имеет функцию глобального поиска, тем самым он является более предпочтительным для идентификации большого числа параметров.
Будем рассматривать задачу в двумерной постановке. Схему расчетной области можно посмотреть на рис. 1, а именно прямоугольную область, отмеченную на схеме красным цветом. Включения будут играть роль адсорбентов, на поверхности которых (Г,), будут происходить гетерогенные реакции, а именно: адсорбция и десорбция. Гомогенные реакции будут происходить в порах (О). Будем считать, что растворы входят на левой границе ) и выходят на правой границе (Гои<). На боковых границах будут стоять условия симметрии
2.1. Гидродинамика. Гидродинамика рассматривается в предположении, что поток влияет на перенос примеси, а примесь на поток не влияет. Вследствие этого скорость просчитывается только один раз и полученный вектор скорости используется для решения нелинейной задачи реагирующего переноса. Поток в пористой среде в масштабе пор считается достаточно медленным и описывается уравнениями Стокса
где и и р — скорость и давление жидкости соответственно, ^ > 0 — вязкость жидкости, которая предполагается постоянной.
2. Постановка задачи
Vp - ^V2u = 0, v ■ u = 0, x e o,
(1) (2)
Рис. 1. Эскиз вычислительной области.
На входной границе задается равномерный поток жидкости:
и п = й, и х п = 0, г £ Гщ. (3)
На выходе задано давление р и условие отсутствия тангенциальных сил:
р — (тп ■ п = р, сгп х п = 0, х € Г0„4, (4)
где а — тензор вязких напряжений, который выглядит следующим образом:
а = ^(Уи + (Уи)т).
На границе включений ставятся условия твердых стенок:
и ■ п = 0, и х п = 0, х £ Гв. (5)
На границах симметрии используются условия идеального скольжения:
и ■ п = 0, ап х п = 0, х £ Г^т. (6)
2.2. Перенос примесей. Буквы А и В будем использовать для обозначения двух растворов. Схемой гомогенной реакции был выбран кубический автокатализ [12], когда
А + 2В ^ 3В, г = каЪаСаС2ь,
где г — скорость реакции. На поверхности адсорбентов (Г8) происходит изотермическая гетерогенная реакция, которая описывается изотермой Ленгмюра
[25]:
А ^ В, Г = ка^СЛ 1 —
— кйе8ТО.
Управляющие параметры: каЪз — коэффициент абсорбции на единицу объема, то > 0 является максимально возможной адсорбированной поверхностной
Г
Г
Г
Г
П
т
т
оо
концентрацией, ка^8 — коэффициент адсорбции на единицу измерения, к^ез -коэффициент десорбции на единицу времени.
Перенос примеси определяется уравнениями конвекции-диффузии
дС
—^ + У-{иСа) = 0(У2Са-каЬзСаС1 х&П, ¿>0, (7)
дЬ
^ + У -{иСь) = ОьУ2Сь + каЬ,СаС2ъ, ¿>0, (8)
со следующими граничными условиями:
Сг = С1, х&Тт, (9)
DlУCl ■ п = 0, г = [а,Ь], х е Т31т и Тот- (10)
Сам конвективный поток через границу выхода разрешен неявно приведенными выше уравнениями. Гетерогенные реакции, которые возникают на поверхности включений Г,, удовлетворяют закону сохранения массы. В данном конкретном случае это означает, что изменение адсорбированной поверхностной концентрации равно потоку от жидкости к поверхности. Массовый баланс на поверхности включений для обоих растворов описывается так:
дт
- = ва\гса-п, жег,, (И)
дт
— = -Въ\ГСъ-п, жег,, (12)
дЬ
где т — концентрация адсорбированного вещества на поверхности, Di > 0 — коэффициент диффузии, С — концентрация веществ А и В соответственно (г = [а,Ь]).
Концентрация адсорбированного вещества на поверхности описывается изотермой Ленгмюра:
дт т
— = каЛаСа 1--- кЛеат, (13)
дЬ \ тж)
где тж является максимально возможной концентрацией на поверхности включений, а ка^8 и клеэ являются параметрами скорости адсорбции и десорбции соответственно.
Начальные условия выглядят следующим образом и считаются известными:
С(х, 0) = 0, г = [а, Ь], х е П, (14)
т(х, 0) = 0, х е Г,. (15)
2.3. Безразмерная форма уравнений. В задачах идентификации параметров является стандартным подход, когда управляющие уравнения рассматриваются в безразмерной форме. Обезразмеривание будет происходить по высоте вычислительной области I для масштабирования пространственных масштабов, скорость будет масштабироваться через входную скорость й, а концентрация будет масштабироваться по входным концентрациям Са и Сь-
Уравнения Стокса (1), (2) и граничные условия (3), (4) принимают следующий вид:
Ур — У2и = 0, (16)
V- и = 0, х е О, (17)
и ■ п =1, и х п = 0, х е Гг„, (18)
р — ап ■ п = 0, ап х п = 0, х е Гои4, (19)
граничные условия (5) и (6) остаются неизменными.
Уравнения (7), (8) и граничные условия (9)-(12) с гетерогенной реакцией (13) принимают следующий вид:
дС 1
+ V • (иСа) = - БааЬ8 СаСь2, (20)
^+У-(иСь) = -^У2Сь+ВиаЬзСаС1 х€П, ¿>0, (21) д£ Реь
Са = 1, Сь = 1, х е Г„, (22)
УСа ■ п = 0, УСь ■ п = 0, X е Г8гт и Гоии (23)
дт
— = УСа-п, жеГ8, (24)
дт
— = -УСь-п, жег,, (25)
дт / ■т- \
— = Ба^ Са - — ^ - Ба^ то, ж е Г,. (26)
Здесь числа Пекле для концентраций обоих растворов:
й1
Рег = —, г = [а, о],
числа Дамколера, которые контролируют кинетику реакции:
,. _ ,. _ ка(13 _ к^ез1
иа аЪв = _ , иа аЛв = _ ; = И ;
и и и
и безразмерный параметр, отвечающий за максимально возможную концентрацию раствора на поверхности адсорбентов:
1с
3. Численное исследование
Для моделирования процессов, проходящих в пористых средах, можно применять различные методы дискретизации по пространству. Метод конечных элементов на текущий момент является эталоном численных методов, объединяя в себе научный и инженерные подходы [26], поскольку хорошо подходит для решения задач на неструктурированной сетке с возможностью измельчения
локального размера ячеек в отдельно взятых подобластях. Единственным недостатком метода конечных элементов можно назвать сложность программной реализации в отличие от других численных методов. На сегодняшний день этим недостатком можно пренебречь, поскольку появились вычислительные библиотеки, которые берут на себя рутинные моменты метода, давая исследователю возможность больше концентрироваться на решаемой задаче [20].
3.1. Вариационная формулировка. Для решения задачи методом конечных элементов нужно переписать уравнения в вариационной формулировке, используя стандартный метод Галеркина. Для скорости определяем пространства пробных и тестовых векторных функций:
V = {u G H: u ■ n =1, u х n = 0 на rin, u = 0 на Ts, u ■ n = 0 на Tsim},
V = {v G Hx(0) : v = 0 на Tm U rs, v ■ n = 0 на Tsim}.
Для давления определяем следующее пространство:
Q = {p G Ь2(Щ, Q = {q G Ь2(П) : q = 0 на Гои<}.
Умножив уравнения на соответствующие тестовые функции и интегрируя по частям, получаем вариационную задачу: найти (u,p) G V х Q, удовлетворяющие граничным условиям (18) и (19), такие, что
У (Vu ■ Vv - V- v p + V - u q) dx = 0 V(v, q) G V х Q. (27)
n
Для нестационарной задачи переноса раствора при гомогенно-гетерогенной реакции определяем следующие пространства:
Sa = {Ca G Hx(0) : Ca = 1 на Гг„}, Sa = {sa G H^fi): sa = 0 на Гг„}, (28)
Sb = {Cb G Hx(0) : Cb = 1 на Гг„}, Sb = {sb G Hx(0) : sb = 0 на Гг„}, (29)
Z = {m G Ь2(Щ, Z = {j G ¿2(0): j = 0 на Гт}. (30)
Необходимо ввести дискретизацию по времени. Будем использовать симметричную схему Кранка — Николсон [27]. Для удобства введем следующее обозначение:
I (ПП+1 I (пп
= --Pn = v(tn), tn =пт, та = 0,1...,
где т — шаг по времени. Найти (Ca,Cb,m) G Sa х Sb х Z, удовлетворяющие граничным условиям (22)—(26):
/Cn+1 _c п f i 1 /* 1
-a-sa d,x + Ca + ~2u .Vsadx + — vcr+i • Vsa d®
n n n
+ i Da absCa + 1 sadx+ j {u ■ n)Ca+* sads
/ТП+1 _ /-1П г 1 1 /* 1
+ / -ь-зь г1х + / С™+2 м • г1х + — / УС™+2 • г1х
- j Ъ&аЬа Са +2 +3Ъ(1х + j (и ■ п)С^+ 2 в^З
О Г0„4
утоп+1_топ у
+ / ---3 аз- Ба ас1з Са 2 I 1--) з аз
г3 Г,
+ J Ба^ ^ ¿в = 0. (31)
Для решения вышеописанной задачи будем использовать вычислительную библиотеку ЕЕпЮБ [20]. Пусть С V, Qh С ф — конечно-элементные пространства на смешанных элементах Тейлора — Худа (Р2 — Р\), а С £а,
С $ь и Zh С Z — конечно-элементные пространства на непрерывных элементах Лагранжа первой степени. Нелинейная задача решается итерационным методом Ньютона.
3.2. Численное решение прямой задачи. Для прямой визуализации микроструктуры пористой среды используют экспериментальный метод компьютерной рентгено-томографии. Такие изображения можно использовать как вычислительную сетку, особенно хорошо в связке с методом конечных объемов (воксель — пиксель). Они помогают хорошо понять морфологию пористых структур, однако для оценки изменчивости, связанной с геометрией и составляющими пористых сред, требуются многочисленные эксперименты для анализа статистических характеристик [10]. Синтетические пористые среды воспроизводят основные усредненные характеристики пористых сред, таких как пористость, поверхность включений и т. д. Они позволяют выполнять многократное моделирование в различных условиях и с множественной параметризацией. Также синтетические пористые среды могут быть легко приведены к любым пространственным масштабам. В данной работе была выбрана идеализированная геометрия с периодическими включениями, чтобы свести к нулю эффекты, возникающие из-за геометрии пор [28], например канальное течение, тупиковые поры и т. д. Исследование влияния таких эффектов при идентификации ключевых параметров реагирующего переноса достойно отдельного систематического исследования и в рамках этой работы проводиться не будет. Создание геометрии исследуемой области и генерация вычислительной сетки было произведено с помощью программы СтвИ.
Поскольку рассматриваемый нами процесс происходит в микромасштабе, экспериментально понять скорости реакции в данный момент является сложной и дорогостоящей задачей. Поэтому часто измеряют данные концентрации растворов на выходе из керна. Эти кривые называются кривыми проскока и они могут быть использованы как вспомогательная информация для определения кинетики реакции. Кривая проскока для математической модели вычисляется
I 0.4
] —*— 4760 узлов 18743 узлов —*— 72745 узлов
-)-1-
сь
Г
1 —*— 4760 узлов 18743 узлов —72745 узлое
-,-^
20 I
(а)
20 I
(Ъ)
Рис. 2. Исследования влияния на кривые проскока: (а) шага по сетке, (Ъ) шага по времени.
следующим образом:
СО
опЬ
/ с(х,г, в)^
Гout
Гout
г =
М],
(32)
где в = {Ре, М, БааЬ5} — вектор параметров. В итоге получаем
средние значения концентраций на границе выхода в каждый момент времени.
На рис. 2(а) и 2(Ь) представлены кривые проскока, которые получены на серии сгущающихся сеток по пространству и по времени соответственно. Видно, что на кривую не сильно влияет количество элементов, поэтому выберем сетку с 4760 узлами. В случае с шагом по времени заметно небольшое отличие, поэтому т был выбран равным 0.2.
Расчетное поле концентраций для обоих веществ в разные временные интервалы можно увидеть на рис. 3. Можно заметить, что с момента времени £ = 8 гомогенная реакция происходит заметнее: концентрация вещества Сь увеличи-
(Ь)
Рис. 3. Расчетное поле концентрации при: (а) £ = 4, (Ь) £ = 8.
вается, поскольку молекулы вещества А превращаются в молекулу вещества В. Результаты скорости и давления идентичны, в силу неучета, предыдущей работе [17].
4. Обратная задача химической кинетики
Представляет интерес решение обратной задачи химической кинетики, которая нацелена на идентификацию определяющих параметров кинетики реакции. В общем, необходимо идентифицировать пять параметров по экспериментальным данным кривых проскока обоих веществ. Существует много подходов решения обратных коэффициентных задач. Выделим основные классы методов, вокруг которых больше всего исследований: градиентные, стохастические, статистические и детерминированные [17]. В этой работе будет рассмотрен относительно молодой по сравнению с остальными методами, но набирающий популярность подход под названием метаэвристика.
4.1. Чувствительность модели к параметрам. Процедура анализа чувствительности модели к параметрам является типичной при решении обратной задачи, поскольку необходимо понять, какие параметры сильнее влияют на кинетику реакции. Из рис. 4(а) видно, что параметр скорости гомогенной реакции существенно влияет на кинетику реакции: чем больше значение, тем быст-
Рис. 4.
(а) м.
Чувствительность модели к параметрам: (а) БааЬв, (Ь) (с)
1.0 /•"»
10.6 I --1
рее происходит переход вещества А в В. На рис. 4(Ь)-4^) показаны ключевые параметры гетерогенной реакции. Параметр Баа^ отвечает за адсорбцию, чем больше значение, тем активнее вещество прилипает к поверхности адсорбента, Ба^ отвечает за десорбцию — процесс, обратный адсорбции, а М отвечает за максимально возможную концентрацию на поверхности адсорбента: чем выше значение, тем дольше адсорбент будет копить депонированную массу перед тем, как начнется процесс десорбции. По рис. 4(с) видно, что параметр, отвечающий за десорбцию, мало влияет на кинетику реакции в выбранном диапазоне параметров. На рис. 5(а),(Ь) показано влияние числа Пекле на кинетику реакции: чем меньше значение, тем больше диффузия реагирующего переноса. Видно, что разные числа Пекле для обоих растворов не оказывают влияния друг на друга.
4.2. Задача идентификации параметров. Будем считать, что имеются экспериментальные данные в виде кривых проскока для каждого вещества, по которым необходимо идентифицировать скорости реакций. Подобные обратные задачи часто рассматривают в виде задачи оптимизации, где необходимо минимизировать некоторый функционал. Будем минимизировать следующий
функционал невязки:
T
J (0) = 1 ((
Coaut - Co)2 + (C0ut - Cb)2) dt,
(33)
0
где Са и Сь — экспериментальные измерения для веществ А и В соответственно, С0и1 и С™1 — кривые проскока для заданных в, которые являются средним значением концентрации веществ А и В на границе выхода и вычисляются по формуле (32). Следует отметить, что данный функционал невязки можно использовать при условии, что значения кривых проскока для обоих веществ одного порядка. В противном случае необходимо использование весов либо предварительно их отмасштабировать на один порядок.
Предположим, что концентрацию каждого вещества будут измерять разные датчики, тем самым они становятся независимыми. Датчики всегда имеют свои погрешности в каких-то заданных диапазонах и они обычно известны. Поскольку есть подходы и методы, которые позволяют убрать шумы в данных, обязательно нужно попробовать идентифицировать параметры на данных без шума. Представляет интерес, когда для экспериментальных данных каждого вещества все-таки присутствует шум, и притом может присутствовать с разной амплитудой для каждого из веществ. Для исследования таких ситуаций будем использовать так называемый синтетический подход генерации экспериментальных данных. Синтетический подход означает, что нужно выбрать из пространства параметров один вектор параметров в, который будет считаться точным решением. Принято, что точное значение для каждого параметра равняется:
Ре = 10.0, Баа^ = 0.005, Ба^ = 0.05, М = 0.05, БааЬ8 = 0.005.
Для простоты исследования и поскольку разные числа Пекле в обоих веществах не оказывают влияния на кинетику друг друга (см. рис. 5(Ь),(с)), берем в расчет частный случай, когда число Пекле взято одинаковым для обоих веществ. Экспериментальные данные генерируются по следующей формуле:
где S — амплитуда шума (представлены данные для 0.01 и 0.05), а ф — случайное число с равномерным распределением в диапазоне [— 1,1]. На рис. 6 представлены описанные синтетические экспериментальные данные с разными амплитудами шума. Далее будем их комбинировать.
4.3. Метаэвристический подход. Метаэвристические методы — относительно новый и быстро развивающийся класс методов оптимизации. Существуют различные метаэвристические алгоритмы, вот несколько наиболее известных: Genetic Algorithms, Particle Swarm Optimization, Differential Evolution и Greedy Randomized Adaptive Search Procedure, подробнее можно ознакомиться
Ci = C(x,t, 0) + Sф, i = [a,b],
(34)
1.0
0,8
— с,
— О. ( '
/
/
-\ /
20
1.0 0.3
I 0.4
— с,
— О.
У
20 е
(а)
(Ъ)
1.0 « 08
£
| 0.6
I 41
зГ 0.4
I
0.2
— С, - Сь
20 t
(с)
Рис. 6. Кривые проскока двух концентраций: (а) без шума, (Ъ) с шумом 8 = 0.01, (с) с шумом 8 = 0.05.
в [29] и по указанным в ней ссылкам. Нами используется Алгоритм Искусственной Пчелиной Колонии (АИПК) [23], который хорошо себя показал в промышленных и инженерных задачах. АИПК имеет функцию глобального поиска, тем самым для идентификации большого числа параметров он является весьма эффективным.
Для идентификации рассмотрим несколько случаев: экспериментальные данные не имеют шума, когда шум в экспериментальных данных одинаков для обоих веществ (6 = 0.01 и 6 = 0.05); когда шум в экспериментальных данных различается в обоих веществах (6а = 0.01 и 6Ь = 0.05 и наоборот). АИПК для всех случаев запускали с одинаковыми управляющими параметрами: число пчел 5, количество итераций 300. Для ознакомления с управляющими параметрами алгоритма можно обратиться к [30]. Критерий останова можно задать в алгоритме, но обычной является практика запуска алгоритма на заранее определенном количестве итераций. Выбор таких параметров для АИПК аргументируется тем, что решаемая система уравнений сама по себе является сложной задачей и если запускать метаэвристический алгоритм (любой мета-эвристический алгоритм, не только АИПК) на слишком большое количество
итераций или с большим количеством пчел, смысл решения задачи может исчезнуть из-за высоких вычислительных затрат. То же самое касается критерия останова. Если задать точность по априорной оценке минимума функционала, алгоритм может произвести слишком большое количество вызовов минимизируемой функции в попытке достичь наилучшего результата, что может повлечь собой большие вычислительные затраты. Следовательно, требуется адекватная оценка имеющихся вычислительных ресурсов.
Для каждого параметра были заданы диапазоны, внутри которых должно быть искомое значение:
Ре = [0,20], Баа^ = [0.0, 0.01], Ба^ = [0.0, 0.1],
М = [0.0,0.1], БааЪв = [0.0,0.01].
При работе с реальными данными такие диапазоны должны быть определены через обоснованное предположение, например, по физическому смыслу. В табл. 1 можно увидеть результаты работы алгоритма: в первом столбце стоят параметры с точными значениями, в остальных указаны идентифицированные параметры для каждого случая, в последней строке показано значение целевого функционала в идентифицированных точках. Видно, что с увеличением шума в экспериментальных данных значение функционала увеличивается, что говорит о том, что он сглаживается. Параметр Ба^ идентифицировался плохо из-за того, что он не сильно влиял на кинетику реакции (см. рис. 4(с)). Слабое влияние этого параметра на характер реагирующего переноса отчасти случилось по причине большого значения параметра М, который ввиду названной причины тоже идентифицировался не очень хорошо. Параметры Баа^ и Бааъ« идентифицировались хорошо на всех случаях. Число Пекле в целом идентифицировалось весьма точно, за исключением случая смешанной амплитуды шума (5а = 0.01 и 5ъ = 0.05). Это можно объяснить тем, что 5 пчел, скорее всего, недостаточно для такого случая. В табл. 2 показаны относительные погрешности для каждого рассмотренного случая в процентах. Сравнительные графики можно увидеть на рис. 7(а), где показаны кривые проскока по идентифицированным параметрам. Графики наложились друг на друга, поэтому нужно обратить внимание на рис. 7(Ь), где показана относительная погрешность в Ь2-норме в каждый момент времени в логарифмической шкале. Погрешность падает до 10-5, особенно в интервалах времени от 10 до 25, когда кривая проскока быстро меняется. Как можно видеть из табл. 2, погрешность для всех случаев получилась меньше 1%, что является вполне приемлемым результатом. Естественно, можно увеличить количество пчел и число итераций для получения параметров с большей точностью.
5. Заключение
В работе реализована модель для описания кинетики гомогенно-гетерогенной реакции и выполнена идентификация параметров скоростей реакции. Вы-
Таблица 1. Идентифицированные параметры для различного уровня шума
Без шума 5 = 0.01 5 = 0.05 5а = 0.01 и 5Ъ = 0.05 5а = 0.05 и 5Ъ = 0.01
Ре = 10.0 О.ЭЭбе I 01 1.045е I 01 1.025е I 01 1.447е I 01 1.095е I 01
Баа^ = 0.005 4.974е-03 4.829е-03 4.927е-03 5.211е-03 4.892е-03
Ба*» = 0.05 9.462е-02 7.722е-02 8.351е-02 8.478е-02 7.429е-02
М = 0.05 9.75е-02 7.490е-02 7.389е-02 7.036е-02 7.751е-02
Бааьв = 0.005 5.053е-03 5.150е-03 4.819е-03 4.942е-03 4.817е-03
3 (0) 3.56е-03 5.053е-02 2.532е-01 1.817е-01 1.817е-01
Таблица 2. Относительная погрешность кривых проскока по найденным параметрам для каждого вещества
Без шума 5 = 0.01 5 = 0.05 5а = 0.01 и 5Ъ = 0.05 5а = 0.05 и 5Ъ = 0.01
¿ае( (%) 0.047 0.049 0.18 0.12 0.19 0.17 0.53 0.43 0.47 0.42
- Са - Сь * Са 6 = 0.01 * Сь 5 = 0.01 • Са 6 = 0.05 • Сь 6 = 0.05 X Са 6 = 0.01 х Сь 6 = 0.05 ▼ Са 6-0.05 ▼ Сь 6 = 0.01 ♦ Са без шума _______
(а) (Ь)
Рис. 7. Результаты идентификации: (а) кривые прорыва двух концентраций, (Ь) относительная погрешность в ¿2 норме в отдельные моменты времени.
числения представлены в двумерном случае на идеализированной геометрии. Полученные результаты формулируются следующим образом.
1. Реализована модель гомогенно-гетерогенной реакции в масштабе пор, основанная на уравнениях Стокса и уравнениях конвекции-диффузии-реакции с граничным условием третьего рода на границах включений. Гомогенная реакция описывается как кубический автокатализ на всем поровом пространстве, а гетерогенная реакция описывается изотермой Ленгмюра.
2. Численное решение задачи производится методом конечных элементов на кусочно-линейных элементах. Дискретизация по времени проводится симметричной схемой Кранка — Николсон. Нелинейная задача решается итерационным методом Ньютона. Массоперенос смоделирован с рассчитанным полем скорости (массоперенос не влияет на скорость). Проведены исследования на
серии сеток для подтверждения сходимости, аналогичные исследования проведены для разных шагов по времени. Дополнительно проведен анализ чувствительности модели к управляющим параметрам реагирующего переноса через пористую среду.
3. Проведена идентификация ключевых параметров, характеризующих реагирующий перенос на основе двух кривых проскока двух разных веществ. Были рассмотрены данные без шума и зашумленные данные с разными амплитудами шума, включая смешанные амплитуды. Для многомерной обратной задачи был применен метаэвристический Алгоритм Искусственной Пчелиной Колонии, который показал хорошую эффективность при достаточно малой вычислительной цене. Параметры были идентифицированы с погрешностью менее 1% для кривых проскока.
ЛИТЕРАТУРА
1. Liu P., Yao J., Couples G. D., Ma J., Iliev O. 3-D modelling and experimental comparison of reactive flow in carbonates under radial flow conditions // Sci. Rep. 2017. V. 7, N 1. P. 1—10.
2. Greiner R., Prill T., Iliev O., van Setten B. A. A. L., Votsmeier M. Tomography based simulation of reactive flow at the micro-scale: Particulate filters with wall integrated catalyst // Chem. Eng. J. 2019. V. 378. Article ID 121919.
3. Myers R. The basics of chemistry. Greenwood Publ. Group, 2003.
4. Yang Q., Zhong Y., Li X.-M., et al. Adsorption-coupled reduction of bromate by Fe(II)-Al(III) layered double hydroxide in fixed-bed column: Experimental and breakthrough curves analysis // J. Ind. Eng. Chem. 2015. V. 28. P. 54-59.
5. Wartha E., Birkelbach F., Bosenhofer M., et al. Enhanced kinetic model identification for gas?solid reactions through computational fluid dynamics // Chem. Eng. J. 2022. V. 430. 132850.
6. Chen L., He A., Zhao J., et al. Pore-scale modeling of complex transport phenomena in porous media // Progr. Energy Combustion Sci. 2022. V. 88. 100968.
7. Sadhukhan S., Gouze P., Dutta T. Porosity and permeability changes in sedimentary rocks induced by injection of reactive fluid: A simulation model //J. Hydrology. 2012. V. 450. P. 134-139.
8. Miller K., Vanorio T., Keehm Y. Evolution of permeability and microstructure of tight carbonates due to numerical simulation of calcite dissolution //J. Geophys. Res., Solid Earth. 2017. V. 122, N 6. P. 4460-4474.
9. Shimura T., Jiao Z., Shikazono N. Evaluation of nickel-yttria stabilized zirconia anode degradation during discharge operation and redox cycles operation by electrochemical calculation // J. Power Sources. 2016. V. 330. P. 149-155.
10. Mosser L., Dubrule O., Blunt M. J. Reconstruction of three-dimensional porous media using generative adversarial neural networks // Phys. Rev. E. 2017. V. 96, N 4. Article ID 043309.
11. Liu Y. Some consideration on the Langmuir isotherm equation // Colloids Surfaces, A, Physic-ochem. Eng. Aspects. 2006. V. 274, N 1-3. P. 34-36.
12. Merkin J. H. A model for isothermal homogeneous-heterogeneous reactions in boundary-layer flow // Math. Comput. Model. 1996. V. 24, N 8. P. 125-136.
13. Kameswaran K., Shaw S., Sibanda P., Murthy P. V. S. N. Homogeneous-heterogeneous reactions in a nanofluid flow due to a porous stretching sheet // Int. J. Heat Mass Transfer. 2013. V. 57, N 2. P. 465-472.
14. Khan M. I., Waqas M., Hayat T., Alsaedi A. A comparative study of casson fluid with homogeneous-heterogeneous reactions //J. Colloid Interface Sci. 2017. V. 498. P. 85-90.
15. Waqas M. A mathematical and computational framework for heat transfer analysis of ferromagnetic non-newtonian liquid subjected to heterogeneous and homogeneous reactions // J. Magnet. Magn. Mater. 2020. V. 493. Article ID 165646.
16. Mozaffari M., Kordzadeh-Kermani V., Ghalandari V., et al. Adsorption isotherm models: A comprehensive and systematic review (2010-2020) // Sci. Total Environ. 2022. V. 812. Article ID 151334.
17. Grigoriev V. V., Iliev O., Vabishchevich P. N. Computational identification of adsorption and desorption parameters for pore scale transport in periodic porous media //J. Comput. Appl. Math. 2020. V. 370. Article ID 112661.
18. Grigoriev V. V., Vabishchevich P. N. Bayesian estimation of adsorption and desorption parameters for pore scale transport // Math., MDPI. 2021. V. 9, N 16. P. 1-16.
19. Grigoriev V. V., Savvin A. V. Numerical study of the influence of the electrokinetic effect on the growth of an oxide film at the pore scale // AIP Conf. Proc. 2022. V. 2528, N 1. Article ID 020047.
20. Aln^s M. S., Blechta J., Hake J., et al. The FEniCS project version 1.5 // Arch. Numer. Software. 2015. V. 3, N 100.
21. Liu Q., Li X., Liu H., et al. Multi-objective metaheuristics for discrete optimization problems: A review of the state-of-the-art // Appl. Soft Comput. 2020. V. 93. 106382.
22. Karimi-Mamaghan M., Mohammadi M., Meyer P., et al. Machine learning at the service of meta-heuristics for solving combinatorial optimization problems: A state-of-the-art // Eur. J. Oper. Res. 2022. V. 296, N 2. P. 393-422.
23. Karaboga D., Akay B. A comparative study of artificial bee colony algorithm // Appl. Math. Comput. 2009. V. 214, N 1. P. 108-132.
24. Grigoriev V. V., Iliev O., Vabishchevich P. N. On parameter identification for reaction-dominated pore-scale reactive transport using modified bee colony algorithm // Algorithms. 2022. V. 15, N 1.
25. Kralchevsky P. A., Danov K. D., Denkov N. D. Chemical physics of colloid systems and interfaces. Handbook of surface and colloid chemistry. CRC Press, 1997.
26. Reddy J. N. Introduction to the finite element method. McGraw-Hill Education, 2019.
27. Samarskii A. A. The theory of difference schemes CRC Press, 2001.
28. Mohring J., Milk R., Ngo A., et al. Uncertainty quantification for porous media flow using multilevel Monte Carlo // Int. Conf. Large-Scale Sci. Comput. Springer, 2015. P. 145-152.
29. Dokeroglu T., Sevinc E., Kucukyilmaz T., Cosar A. A survey on new generation metaheuristic algorithms // Comput. Ind. Eng. 2019. V. 137. Article ID 106040.
30. Song X., Zhao M., Yan Q., Xing S. A high-efficiency adaptive artificial bee colony algorithm using two strategies for continuous optimization // Swarm Evolut. Comput. 2019. V. 50. Article ID 100549.
Поступила в редакцию 3 марта 2023 г. После доработки 29 апреля 2023 г. Принята к публикации 29 мая 2023 г.
Григорьев Василий Васильевич
Лаборатория «Вычислительные технологии моделирования многофизичных и многомасштабных процессов криолитозоны». Северо-Восточный федеральный университет имени М. К. Аммосова, ул. Кулаковского 42, Якутск 677000;
Северо-Кавказский центр математических исследований, ул. Пушкина, 1, Ставрополь 355000 [email protected]
Математические заметки СВФУ Апрель—июнь, 2023. Том 30, № 2
UDC 519.63
IDENTIFICATION OF HOMOGENEOUS-HETEROGENEOUS PORE-SCALE REACTION RATES IN POROUS MEDIA V. V. Grigoriev
Abstract: This paper presents a model of homogeneous-heterogeneous reaction in the pore scale based on Stokes equations, convection-diffusion-reaction equations with the Robin boundary condition at the inclusion boundaries. The homogeneous reaction is described as cubic autocatalysis on the whole pore space, and the kinetics of the heterogeneous reaction is described by the Langmuir isotherm. Numerical solution of the problem is carried out by the finite element method on piecewise linear elements. The Crank-Nicholson scheme is used for discretization in time. The nonlinear problem is solved using Newton's iteration method. The mass transfer is simulated with a calculated velocity field. In addition, a sensitivity analysis of the model to the parameters has been carried out to study their influence on the reactive transport through the porous medium. A numerical solution for the inverse problem, namely, identification of key parameters characterizing the reactive transport based on two breakthrough curves of two different solutions is presented. Noisy measurements with different noise amplitudes including mixed amplitudes were considered. For approximate solution of the multidimensional inverse problem the metaheuristic Artificial Bee Colony Algorithm was applied and showed good efficiency at rather low computational cost.
DOI: 10.25587/SVFU.2023.74.45.008
Keywords: homogeneous-heterogeneous reaction, porous media, pore scale, parameter identification, finite element method.
REFERENCES
1. LiuP., YaoJ., Couples G. D., Ma J., and Iliev O., "3-D modelling and experimental comparison of reactive flow in carbonates under radial flow conditions," Sci. Rep., 7, 1—10 (2017).
2. Greiner R., Prill T., Iliev O., van Setten B. A. A. L., and Votsmeier M., "Tomography based simulation of reactive flow at the micro-scale: Particulate filters with wall integrated catalyst," Chem. Eng. J., 378, article No. 121919 (2019).
3. Myers R., The Basics of Chemistry, Greenwood Publ. Group (2003).
4. Yang Q., Zhong Y., Li X.-M., et al., "Adsorption-coupled reduction of bromate by Fe(II)-Al(III) layered double hydroxide in fixed-bed column: Experimental and breakthrough curves analysis," J. Ind. Eng. Chem., 28, 54-59 (2015).
5. Wartha E., Birkelbach F., Bösenhofer M., et al., "Enhanced kinetic model identification for gas-solid reactions through computational fluid dynamics," Chem. Eng. J., 430, article No. 132850 (2022).
6. Chen L., He A., Zhao J., et al., "Pore-scale modeling of complex transport phenomena in porous media," Progress Energy Combust. Sci., 88, article No. 100968 (2022).
7. Sadhukhan S., Gouze P., and Dutta T., "Porosity and permeability changes in sedimentary rocks induced by injection of reactive fluid: A simulation model," J. Hydrol., 450, 134-139 (2012).
© 2023 V. V. Grigoriev
8. Miller K., Vanorio T., and Keehm Y., "Evolution of permeability and microstructure of tight carbonates due to numerical simulation of calcite dissolution," J. Geophys. Res., Solid Earth, 122, No. 6, 4460-4474 (2017).
9. Shimura T., Jiao Z., and Shikazono N., "Evaluation of nickel-yttria stabilized zirconia anode degradation during discharge operation and redox cycles operation by electrochemical calculation," J. Power Sources, 330, 149-155 (2016).
10. Mosser L., Dubrule O., and Blunt M. J., "Reconstruction of three-dimensional porous media using generative adversarial neural networks," Phys. Rev. E, 96, No. 4, article No. 043309 (2017).
11. Liu Y., "Some consideration on the Langmuir isotherm equation," Colloids Surfaces, A, Physic-ochem. Eng. Aspects, 274, No. 1-3, 34-36 (2006).
12. Merkin J. H., "A model for isothermal homogeneous-heterogeneous reactions in boundary-layer flow," Math. Comput. Model., 24, No. 8, 125-136 (1996).
13. Kameswaran K., Shaw S., Sibanda P., and Murthy P. V. S. N., "Homogeneous-heterogeneous reactions in a nanofluid flow due to a porous stretching sheet," Int. J. Heat Mass Transfer, 57, No. 2, 465-472 (2013).
14. Khan M. I., Waqas M., Hayat T., and Alsaedi A., "A comparative study of casson fluid with homogeneous-heterogeneous reactions," J. Colloid Interface Sci., 498, 85-90 (2017).
15. Waqas M., "A mathematical and computational framework for heat transfer analysis of ferromagnetic non-newtonian liquid subjected to heterogeneous and homogeneous reactions," J. Magnet. Magn. Mater., 493, article No. 165646 (2020).
16. Mozaffari M., Kordzadeh-Kermani V., Ghalandari V., et al., "Adsorption isotherm models: A comprehensive and systematic review (2010-2020)," Sci. Total Environment, 812, article No. 151334 (2022).
17. Grigoriev V. V., Iliev O., and Vabishchevich P. N., "Computational identification of adsorption and desorption parameters for pore scale transport in periodic porous media," J. Comput. Appl. Math., 370, article No. 112661 (2020).
18. Grigoriev V. V. and Vabishchevich P. N., "Bayesian estimation of adsorption and desorption parameters for pore scale transport," Math., MDPI, 9, No. 16, 1-16 (2021).
19. Grigoriev V. V. and Savvin A. V., "Numerical study of the influence of the electrokinetic effect on the growth of an oxide film at the pore scale," AIP Conf. Proc., 2528, No 1, article No. 020047 (2022).
20. Aln^s M. S., Blechta J., Hake J., et al., The FEniCS Project Version 1.5, Arch. Numer. Software, 3, No. 100 (2015).
21. Liu Q., Li X., Liu H., et al., "Multi-objective metaheuristics for discrete optimization problems: A review of the state-of-the-art," Appl. Soft Comput., 93, article No. 106382 (2020).
22. Karimi-Mamaghan M., Mohammadi M., Meyer P., et al., "Machine learning at the service of meta-heuristics for solving combinatorial optimization problems: A state-of-the-art," Eur. J. Oper. Res., 296, No. 2, 393-422, (2022).
23. Karaboga D. and Akay B., "A comparative study of artificial bee colony algorithm," Appl. Math. Comput., 214, No. 1, 108-132 (2009).
24. Grigoriev V. V., Iliev O., and Vabishchevich P. N., "On parameter identification for reaction-dominated pore-scale reactive transport using modified bee colony algorithm," Algorithms, 15, No. 1 (2022).
25. Kralchevsky P. A., Danov K. D., and Denkov N. D., "Chemical Physics of Colloid Systems and Interfaces," in: Handbook of Surface and Colloid Chemistry, CRC Press (1997).
26. Reddy J. N., Introduction to the Finite Element Method, McGraw-Hill Education (2019).
27. Samarskii A. A., The Theory of Difference Schemes, vol. 240, CRC Press (2001).
28. Mohring J., Milk R., Ngo A., et al., "Uncertainty quantification for porous media flow using multilevel Monte Carlo," in: Int. Conf. Large-Scale Scientific Computing, pp. 145-152, Springer (2015).
29. Dokeroglu T., Sevinc E., Kucukyilmaz T., and Cosar A., "A survey on new generation meta-heuristic algorithms," Comput. Ind. Eng., 137, article No. 106040 (2019).
30. Song X., Zhao M., Yan Q., and Xing S., "A high-efficiency adaptive artificial bee colony algorithm using two strategies for continuous optimization," Swarm Evolut. Comput., 50,
article No. 100549 (2019).
Submitted March 3, 2023 Revised April 29, 20023 Accepted May 29, 20023
Vasiliy V. Grigoriev
Laboratory of Computational Technologies
for Modeling Multiphysical and Multiscale Permafrost Processes, Ammosov North-Eastern Federal University, 42 Kulakovskii Street, Yakutsk, 677000 Russia; North Caucasus Centre for Mathematical Research, 1 Pushkin Street, Stavropol, 355000 Russia [email protected]
