ИДЕНТИФИКАЦИЯ ПАРАМЕТРОВ ТРАНСПОРТА ПРИМЕСИ ЧЕРЕЗ ПОРИСТУЮ СРЕДУ С УЧЁТОМ ЗАКУПОРКИ Текст научной статьи по специальности «Математика»

ВЕСТНИК ПЕРМСКОГО УНИВЕРСИТЕТА

2021

• ФИЗИКА •

Вып. 3

УДК 532.546; 532.5.013.4 PACS 47.56.+Г, 47.55.pd

Идентификация параметров транспорта примеси через пористую среду с учётом закупорки

М. Р. Хабин1'2^, Б. С. Марышев1'2!

1 Пермский государственный национальный исследовательский университет, Пермь, Россия

2 Институт механики сплошных сред УрО РАН, Пермь, Россия

t email: mikhail.khabin@psu.ru

* email: bmaryshev@mail.ru

Статья посвящена определению параметров транспорта примеси в пористой среде, гидродинамическое сопротивление которой меняется вследствие оседания частиц примеси на стенках пор (закупорки). В большинстве случаев экспериментальное измерение параметров такого транспорта не является возможным. При нахождении таких параметров для исследуемой модели транспорта может быть решена обратная задача. Обратная задача решается для некоторого набора экспериментальных данных, которые можно соотнести с результатами, полученными при моделировании, с помощью BFGS-метода в сочетании с методом присоединённой функции. В данной работе моделируется прокачка фиксированного объёма примеси через рабочую область, заполненную пористой средой, при постоянном перепаде давления на её концах. В качестве «экспериментального» набора данных берутся расход и концентрация примеси на выходе из рабочей области, полученные для исследуемой модели с заданными параметрами. Основным из подходов для описания такого транспорта является MIM (mobile/immobile media) подход, заключающийся в разделении общей концентрации примеси на концентрацию мобильную и немобильную. Для небольших концентраций примеси может быть применена нелинейная MIM модель. Её кинетическое уравнение содержит концентрацию насыщения немобильной компоненты, при достижении которой прекращается процесс адсорбции частиц. Закупорка в модели учитывается изменением проницаемости пористой среды вследствие изменения пористости, которая, в свою очередь линейно, уменьшается с ростом концентрации немобильной компоненты. Для рассматриваемой модели была решена обратная задача по нахождению заданного набора параметров. Приведены результаты тестовой работы алгоритма решения обратной задачи. Показано, что алгоритм сходится к заданным параметрам за небольшое число итераций. Предполагается, что данный алгоритм будет использоваться для обработки данных эксперимента с цилиндрической рабочей областью, заполненной стеклянным гранулятом в качестве пористой среды.

Ключевые слова: пористая среда; обратная задача; закупорка

Поступила в редакцию 15.08.2021; после рецензии 31.08.2021; принята к опубликованию 31.08.2021

Identification of the parameters of transport through porous media with clogging

M. R. Khabin12t, B. S. Maryshev12i

1 Perm State University, Perm, Russia

2 Institute of Continuous Media Mechanics UB RAS, Perm, Russia t email: mikhail.khabin@psu.ru

* email: bmaryshev@mail.ru

© Хабин М. Р., Марышев Б. С., 2021

распространяется на условиях лицензии

Creative Commons Attribution 4.0 International (CC BY 4.0).

The article is devoted to the determination of the parameters of impurity transport in a porous medium, the hydrodynamic resistance of which changes due to the settling of impurity particles on the pore walls (blockages). In most cases, experimental measurement of the parameters of such transport is not possible. When finding such parameters for the transport model under study, the inverse problem can be solved. The inverse problem is solved, for a certain set of experimental data, which can be correlated with the results obtained in modeling, using the BFGS method in combination with the adjoint function method. This paper simulates the pumping of a fixed volume of impurities through a working area filled with a porous medium at a constant pressure drop at its ends. As an "experimental" data set, the flow rate and concentration of impurities at the exit from the working area, obtained for the investigated model with the given parameters, are taken. The main approach for describing such transport is the MIM (mobile / immobile media) approach, which consists in dividing the total concentration of an impurity into a concentration of mobile and nonmobile. For small impurity concentrations, a nonlinear MIM model can be applied, the kinetic equation of which contains the saturation concentration of the immobile component, upon reaching which the adsorption of particles stops. Blockage in the model is taken into account by a change in the permeability of a porous medium due to a change in its porosity, which in turn decreases linearly with an increase in the concentration of the immobile component. For the model under consideration, the inverse problem of finding a given set of parameters was solved. The results of the test operation of the algorithm for solving the inverse problem are presented. It is shown that the algorithm converges to the specified parameters in a small number of iterations. It is assumed that this algorithm will be used to process the experimental data with a cylindrical working area filled with glass granulate as a porous medium.

Keywords: porous media; inverse problem; blockage

Received 15.08.2021; revised 31.08.2021; accepted 31.08.2021

doi: 10.17072/1994-3598-2021-3-44-55 1. Введение

1.1. Актуальность работы

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

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

Модель, параметры которой определяются в данной работе, описывает стандартную постановку эксперимента, [1,2] по измерению транспортных параметров пористой среды. Однако помимо изучения влияния структуры среды на концентрационный профиль (адсорбции примеси), учитывается и обратное влияние оседания примеси на фильтрационные свойства - закупорка. Последний

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

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

Существует множество методов решения оптимизационной задачи. Стандартная классификация таких методов исходит из количества членов, оставляемых в ряде Тейлора при разложении целевой функции около точки минимума по параметрам модели. Так, методы, основанные на расчете целевой функции в точке без расчета её производных, называются методами нулевого порядка или методами прямого поиска. Это такие методы, как метод Хука-Дживса, метод касательных, метод деформируемого многогранника [4]. Методы, для реализации которых требуется вычисление первых производных по параметрам модели, называются методами первого порядка или методами спуска, к ним относится метод наискорейшего спуска [5] и метод сопряженных градиентов [6]. Квази-Ньютоновские методы второго порядка основываются на информации о самой функции, ее первой и второй производных. Наиболее известные методы такого класса - это метод Девидона - Флетчера - Пауэлла (DFP) [7] и Бройдена - Флетчера - Гольдфарба - Шанно (BFGS) [8]. Последний класс методов является наиболее точным и быстросходящимся, однако существует проблема при вычислении градиента целевой функции в пространстве параметров модели. Ошибки при их вычислении сильно замедляют сходимость метода. Наиболее точным способом вычисления таких градиентов является метод присоединенной функции [9,10]. Настоящая работа посвящена решению обратной задачи по иден-

тификации параметров переноса примеси фильтрационным потоком в пористой среде. Представляется, что наиболее оптимальным методом её решения является BFGS-метод в сочетании с методом присоединенной функции.

Диффузию в пористой среде обычно описывают с помощью стандартного закона Фика в предположении того, что поток частиц линейно зависит от градиента концентрации. Часто при таком описании возникают расхождения с экспериментальными данными, вследствие чего требуется модификация стандартного подхода. Одним из вариантов такой модификации является MIM (mobile-immobile media) подход [11], в рамках которого предполагается, что общая концентрация примеси может быть разделена на две составляющих: неподвижной фаза и подвижная фаза. Первая связана с частицами, осевшими (или иммобилизованными) на твердый скелет пористой среды, а вторая с примесью, переносимой как диффузионным, так и адвективным потоками. При таком описании требуется дополнительное уравнение, которое называют «кинетическим уравнением» -оно описывает переход частиц между двумя фазами. В зависимости от вида кинетического уравнения, существуют различные модели для описания диффузии в пористой среде.

Стандартная модель MIM [12] - это линейная модель, в которой поток примеси в неподвижную фазу линейно возрастает с ростом концентрации примеси в подвижной фазе и линейно уменьшается с ростом концентрации в неподвижной фазе.

Фрактальная (дробная) модель MIM [13] предполагает, что частицы примеси в пористой среде остаются связанными на протяжении случайных промежутков времени и распределение такой случайной величины не имеет среднего значения (наблюдается субдиффузионный перенос). В этом случае кинетическое уравнение представляет собой линейную зависимость между потоком частиц в неподвижную фазу и дробной производной Ка-путо от концентрации подвижных частиц по времени. Эта модель подтверждена экспериментальными данными, прежде всего, для малых концентраций примеси [14].

Нелинейная MIM модель с насыщением [15] -используется для описания транспорта в пористых средах в присутствии значительного количества растворенного вещества. Предполагается, что концентрация примеси в неподвижной фазе не должна превышать некоторое предельное значение. Этот эффект был рассмотрен в рамках кинетической модели второго порядка, в которой скорость адсорбции линейно зависит от разницы между предельной концентрацией насыщения и концентрацией частиц в неподвижной фазе.

В исследуемых экспериментальных данных наблюдается снижение расхода, предположительно связанное с явлением закупорки пористой сре-

Рис. 1. Разделение выделенного объёма V пористой среды на участки с объёмами: Уц - объём стенок пористой среды, VI - объём жидкости, Vm - объём мобильной компоненты примеси, V - объём немобильной компоненты примеси, Vo - объём порового пространства, свободного от примеси

ды. Как показано в работах [16, 17], закупорка может наблюдаться по причине уменьшения порового пространства за счет оседания (адсорбции) частиц примеси на стенках пор. Известно, что если закупорка происходит достаточно медленно (т.е. не происходит механической закупорки пор, связанной с «затыканием» поры достаточно большой частицей или агломератом таких частиц), то уменьшение порового пространства приводит к снижению проницаемости среды. Зависимость проницаемости среды от объема порового пространства (или пористости) в наиболее общем виде может быть описана в рамках модели Козени-Кармана [18,19]. Закупорка становится заметной только при достаточно больших значениях концентрации примеси в немобильной фазе, что подразумевает использование нелинейной MIM модели с учетом насыщения неподвижной фазы [15]. Задача состоит в том, чтобы для переноса примеси течением в пористой среде определить параметры нелинейной MIM модели на основе экспериментальных данных.

2. Модель транспорта в пористой среде

Для описания транспорта примеси в пористой среде, рассмотрим некоторый объем пористой среды V. Внутри выделенного объема есть пространство, не занятое твердым веществом среды, объём Vo такого пространства называют поровым. Когда через пористую среду протекает смесь, по-ровое пространство заполняется несущей жидкостью и примесью. Предположим, что часть примеси может оседать на стенках скелета пористой среды, тогда объем поры занят объёмами трех компонент: объёмом несущей жидкости V/, объёмом занимаемый примесью осевшей на стенки пор

Vi и объёмом свободной примеси переносимой общим потоком несущей жидкости Vm (рис. 1)

Следовательно, Vo=V/+Vm+Vi, тогда введем новую величину Vp= Vi+Vm - текущий объем поры, т.е. объём пространства, в котором несущая жид-

кость и свободная примесь могут перемещаться. Поделим полученное выражение V0=Vp+Vi на полный объем среды V:

V V V — = —Р- | i

V~ V V'

(1)

Величина V0/V=ф0 является пористостью «чистой» среды без примеси. Отношение Vp/V=ф, в свою очередь, отражает текущую пористость среды, а слагаемое Vi/V=q определяет объемную концентрацию осевших (неподвижных или адсорбированных) частиц примеси. Выразим в новых терминах текущую пористость, тогда равенство (1) перепишется следующим образом:

Ф = Ф0 - q

(2)

Из выражения Vp=Vl+Vm определим объемную концентрацию подвижных частиц с, для это разделим Vp и учтём, что, Vm/Vp=c, тогда

1 = VI + с,

Vp

(3)

где V/Vp является по смыслу объемной концентрацией несущей жидкости в текущем поровом пространстве. Запишем закон сохранения массы примеси.

V^+VI

dt { V V

dt

Ьс + q ) = -div ((J,.). (4)

Поток массы Jс может быть связан только с подвижной составляющей примеси и выражается согласно закону Фика как , где Б - эффективный коэффициент диффузии J с = -ПЧс + ус, у -скорость жидкости в поре или поровая скорость. Подставив поток в (4), получим

d

—(фс + q) = div (фБУс -фус).

(5)

Аналогичным образом можно записать закон сохранения массы жидкости:

f (V)--И V)-v (V1

исходя из (3), тогда (6) перепишется в виде

(6)

5(ф(1 - с))

dt

- div (DV(1 - с) - v (1 - с)) , (7)

div^v) - divu - 0,

(8)

Л

-Vp,

(10)

:(ф( -y

(1 -i

(11)

dqq = a (- q) с - bс,

(12)

суммируя выражения (5) и (7) с учетом (2), а именно ф+q=ф0=const, получим соотношение для фильтрации несжимаемой жидкости:

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

д

— (фс + Ч) = ^ (фС>Ус - ис), Шу и = 0.

Система уравнений (9) незамкнута, и требуется некоторое соотношение, позволяющее связать поле скорости фильтрации с внешними задаваемыми параметрами, а также уравнение, описывающее переход примеси из подвижного состояния в неподвижное, оно должно определять концентрацию Ч. Для связи скорости фильтрации с приложенным внешним давлением воспользуемся законом Дар-си [20]

и = КФФ*

где a, b - коэффициенты адсорбции и десорбции соответственно. Кинетическое уравнение (12) соответствует нелинейной MIM модели [15] и представляет собой эволюционное уравнение для переменной q. Его правая часть имеет слагаемое, содержащее максимально возможное значение концентрации примеси в неподвижной фазе, так называемую концентрацию насыщения немобильной компоненты примеси qo. При достижении равенства q=qo прекращается процесс адсорбции примеси.

Все вышеперечисленные уравнения составляют полную замкнутую систему уравнений, которая описывает транспорт примеси в пористой среде с учетом закупорки:

d

-- (фс + q)- DVфVс + фDАс - uVс,

dq ~dt

u

- a (q0 - q)с - Ьс,

'\Ф)

Л

div u - 0,

Vp,

(13)

к

(Ф)-У

(1 -i

где к(ф) - проницаемость среды, п - коэффициент динамической вязкости. Предположим, что оседание частиц на стенки пор происходит достаточно медленно без существенного изменения формы зерна пористой среды, в этом случае проницаемость может быть описана как однозначная функция пористости. Наиболее популярный способ такого описания дается уравнением Козени-Кармана [18,19]:

ф-ф0- q.

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

d ~dt

dq "dt

div

с + -

\ ф0 J

- DАс +

к(ф) Vp ^с

- a (q0 - q) с - Ьс, ~ (ф), "

-Vp

(14)

к

(ф) -У

- 0,

■ - q)3

(1 -ф0 + q)

где у - параметр Козени-Кармана, зависящий от формы и распределения зерен среды. Переход примеси из подвижного состояния (фазы) в неподвижное будем описывать в рамках MIM подхода [11, 12], что подразумевает разделение примеси на подвижную и неподвижную фазы. В достаточно общей формулировке переход примеси между фазами может быть описан кинетическим уравнением с Ленгмюровской изотермой, характеризующей состояние равновесия [21], а именно:

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

к(ф) ^ к(ф0)

>)=

-к(ф0 )Г(ф )g,

(3 -Фр)

(15)

ф0 (1 -ф0 )

q

с + —

. ф0 J

7

- q) др дс Б д2 с

+ Б я..2'

Т7Ф0 (1 -ф0 + q) дг дг дг

— = а (qo- q )с - Ьс

д2 р-_(3-ф0-?) дР дг2 (ф0 - q)(1 - ф0 - q) дг дг

(16)

= 0,

Рг=0 = Р,

Рг = X = P2, = 7 (фь-^))

•Ур,

Очевидно, что Г(ф0)^ж при ф0^0,1, т.е. когда приближение пористой среды становится неприменимым. Функция Г(ф0) имеет один минимум в интервале ф0 £ (0,1), при ф0 = 3 - -ч/б « 0.55 , при этом Г(ф0*) ~ 10, т.е. Г(ф0)»1. В этом случае необходим учет зависимости проницаемости от концентрации примеси в немобильной фазе, поскольку даже малые вариации пористости приводят к значительным вариациям проницаемости.

3. Постановка прямой задачи

Исследуется одномерное течение двухкомпо-нентной смеси через массив пористой среды. Фильтрация жидкости описывается в приближении Дарси [20], транспорт примеси - с помощью нелинейной М1М модели [15], зависимость пористости от концентрации осевших частиц предполагается линейной. Закупорка описывается в рамках модели Козени-Кармана [18, 19], как взаимно однозначное соответствие между пористостью и проницаемостью среды. Принципиальная схема решаемой задачи представлена на рис. 2.

Л(1 -ф0 - q)

,с1=0 -(ф0 - С0 )Б

1 а )=Г'

10, г > г*

= »1=0 / (г )

(17)

Б * дг

= 0,

Рис. 2. Принципиальная схема одномерного транспорта примеси в рабочей области

Моделируется одномерный транспорт через ограниченную область пористой среды длиной X. Перепад давления между границами области предполагается постоянным (др=Р2-Р1). На левой границе области (входе) в течение интервала времени г £ [0, г*] задается поток примеси с постоянной концентрацией С0, на правой границе задается условие свободного вытекания. Математически задача может быть сформулирована следующим образом:

с (г, г = 0) = 0, q (г, г = 0 ) = 0,

где с, q - объёмные концентрации примеси мобильной компоненты и немобильной компоненты соответственно, q0 - концентрация насыщения немобильной компоненты, р - давление, а, Ь - параметр адсорбции и десорбции, Б - эффективный коэффициент диффузии, ф0 - пористость незагрязнённой среды, и - скорость фильтрации, у - параметр Козени-Кармана, п - динамическая вязкость несущей жидкости.

4. Метод решения прямой задачи

Для решения системы (16), (17) применяется метод конечных разностей второго порядка точности по пространству и первого порядка точности по времени. Решение ищется в области г £ [0, X], г £ [0, Т] на сетке с шагом по пространству Н=Ь/Ы и по времени т=Т/М, следовательно, р(г,гк), с(г„гк), ((г^^ = р^, сгк Х{=Ш, i=1..N, гк=кт, к=1..М. Система уравнений (16), (17), записанная в конечно-разностной форме с использованием неявной схемы [22]:

(3 -ф0- ^) - рк+1 - рк-1 + ( - qk )(1 -ф0 - qk) 2И 2И + рк+1 - 2 рк + рк_1 = 0

^к+1 - qk

к+1 к с, - с.

= а( -qk+1 )ск+1 -

(18)

^ + а( -^+1 )+1 -Ь(^к+1 = ( - ^ )3

рк - рк ск - ск И,+1 У,-1 0,+1 0,-1

Лф0 (1 -фo¡ - ^ ) 2И

2И

к г% к . к

+б5+1 - 2с + с,-1

с граничными и начальными условиями в форме:

и2

Запишем граничные и начальные условия для системы (18):

г=0

т

Ро = Р, А = Рг,

ик ^ХЛ^-С( -р),

Ф (1 -Фо - Ч)

С00(к) Б

(19)

Со ик-1И + Б ' ик-1И + Б ""

в(к ) = |^0к Ч к *, 4 = ^ -1, 10, к > к

с,0 = 0, Ч,0 = 0,

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

е=%-к*.

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

л =

(-ф0 + ^) - & " (1 -Ф,+чк)(- чк) 4 : (л,к +1) - 2Рк +(1 - Л, ) = 0,

в =ут (ф0 -)

Р,к+1 - Р,к-1

, ФФ (1 -ф0 + )2 М2

(-Б + В)с«4!1 + Ффт(«0 -'к))<к" +

(20)

+1- Б- - вк\ с- =Фщк + с

к I ^к+1 — и___к , „к

'0

Я (р ){и} = 0,

(21)

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

ёЕ дЕ дЕ ди

— = — +--. (22)

ар др ди др

Возьмем градиент в пространстве параметров от уравнения (21) и получим

дЯ {и}+я

др [ др

= 0.

(23)

Введём некоторую функцию ¥ и умножим на неё выражение (23)

у * дЯ {и } + У * Я 1 = 0,

др [др

(24)

к+1 = — ск+' + Ч, Ч 1 + Ь- + а—к+1

Полученная система решается в несколько этапов, сначала решается первое и второе уравнения,

к к+1

из них ргк и с, находятся методом прогонки, так как их матрицы имеют трехдиагональный вид. После чего в третье уравнение системы (20) подставляется с,к+1, что позволяет найти Ч,к+1.

5. Обратная задача. Метод присоединенной функции

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

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

у * я {^1 = ^ * Я-1 {У}, (25)

[др ] др 1 '' ' '

где Я-1- оператор обратный Я. Прибавим к правой части равенства (22) левую часть выражения (24) с учетом (25), в результате будем иметь

ёЕ дЕ дЕ ди 1Т, _дЯ,ттЛ — = — +--+ У *—{и} +

ёр др ди др др

ди. др ди. др

■Я-1 {У}= —+ У * 1 ' др

дЯ др

{и}+ (26)

Я-

1 1 ди

Для нахождения функции ¥ положим последнее слагаемое в (26) равным нулю:

Я-1 {У} + — = 0.

ди

(27)

где Я(р) - некоторый оператор, зависящий от вектора параметров р и действующий на наблюдаемую и. Обратная задача формулируется так: «требуется найти такой вектор параметров р, чтобы значение целевой функции Е(и, р) было минимально». В отличие от прямой задачи, которая требует по известному вектору параметров р

Уравнение (27) задаёт собой задачу по нахождению присоединённой функции ¥, сопряженную прямой задаче (21). С учетом решения задачи (27) выражение (26) позволяет рассчитать компоненты градиента целевой функции в виде

^^У^И. (28)

ар др др

6. Метод решения обратной задачи

Систему (22) можно записать в более компактной форме

Ркрк = гк,

] i >

0к к+1 /к . к иск = Х + с-

к+1 к Чг = Я/ ,

(29)

где г1к,/к,як с учётом граничных условий (19) запишутся следующим образом:

-(Ак +1) р, «=1

гк =|0,1 = 2..N - 2 ,

(А»-1 -1),I = N -1

Т +Г ^ - вк 1 ^+^ , г = 1

/к =

икИ + Б

(30)

Ьт к ■ 2 ы 1

—Ч1 ,! = 2..N -1

aтq0 (+1 + С0Ищв^к +1)) + (к ((И + Б) (1 + тЬ )((И + Б) + ат ((с*+1 + С0 Ии0в(к +1))'

„к = атд0ск+1 + дк

о, .

1 + тЬ + атск

,1 = 1..N.

Зададим целевую функцию как сумму квадратов невязок в следующем виде:

е=1:

о» - ок

Qí

К2 ( ск" - Ск" V

"I

ск

^—' Л?

(31)

где 0е»к=0е(Ь,кт) - экспериментально измеренный расход смеси на выходе из рабочей области в момент времени к'т, Се»к=Се(Х,к"т) - значения концентрации примеси на выходе из рабочей области в момент времени к"т, полученные в эксперименте. Наборы натуральных чисел {к'} и {к"}, состоящие из К1 и К2 элементов, соответственно, описывают набор номеров узлов временной сетки, в которых имеются экспериментальные данные, с»к' - концентрация примеси на выходе из рабочей области, полученная в результате расчёта, 0ык - массовый расход смеси на выходе из рабочей области, полученный в результате расчёта. Он может быть вычислен по формуле

о» = (с»р +(1 - с» )р),

)(( - р» -1)

(32)

и» =-

Л

И

единенных функций (ф, щ, х) может быть записана в следующем виде:

Ркук =

,ВТ i

др:

кк- Ук+1 . к Л1

к к V

д/к „к+1

др;

4V+1 .

дрк

дЕ 'дрк

X =

дс

■к к V

„к _ д/ к+1 , д&1 ..к+1 , 1 ...к

дгк

+—V дчГ' Сч^1

к

„к" х

(33)

дЕ дРк к к -т-Т-^ркМ

дО;

дЧ.

К+У+1,

^ -1„к ^ X + к + 1 дЕ

V =~ГГX + М - —.

дск

дск

Последовательно решая систему (33) тем же численным методом, что и задачу (20) получим значения присоединённых функций. После чего на основе выражения (28) найдем интересующие нас компоненты градиента целевой функции в пространстве параметров задачи:

р = {{ а Ь, ч0}

дЕ

др

м

-м

др

др

сМ

(34)

д/к-1 к -—— V +-др

„к-1

др

к к

х, +-с- ^к др

Вычисление компонент градиента целевой функции в пространстве параметров задачи позволяет применять методы второго порядка, как было замечено выше. В нашем случае для нахождения минимума целевой функции используется известный алгоритм, основанный на методе ньютона BFGS (алгоритм Бройдена - Флетчера - Голь-дфарба - Шанно [8]). Результат тестирования его работы для некоторого эталонного набора данных, которые заменяют реальные результаты эксперимента, представлены ниже.

7. Результаты

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

Задача о нахождении компонент градиента целевой приводит к задаче для системы присоединённых функций (ф, щ, х), соответствующих каждой переменной прямой задачи (с,р,ч). Для к=1.М задача решается в обратном (по сравнению с прямой задачей) направлении по времени. В общем виде система разностных уравнений для присо-

Рис. 3. График зависимости целевой функции в логарифмическом масштабе от номера итерации

к

к =1

Рис. 4. Эволюция решения в ходе работы оптимизационного алгоритма. Набор вычисленных данных (пунктирная линия) для концентрации на выходе из рабочей области. Наборы параметров, соответствующих приведенным решениям представлены в таблице. Они сходятся к набору параметров для эталонных данных (сплошная линия), соответствующему а=220, Ь=60,

Б = 1.0■ 10-4, Чо=0.20

Т,с

Т,с

Т,с

Т,с

Рис. 5. Эволюция решения в ходе работы оптимизационного алгоритма. Набор вычисленных данных (пунктирная линия) для потока на выходе из рабочей области. Наборы параметров, соответствующих приведенным решениям представлены в таблице. Они сходятся к набору параметров для эталонных данных (сплошная линия), соответствующему а=220, Ь=60, Б=1.0■ 10-4, Чо=0.20

Полученные значения параметров в зависимости от номера итерации

n=1 n=2 n=10 n=40

a 212 213 213 219.3

b 80 77 77 59.9

D 1.0-10-3 1.0-10-3 1.0-10-3 1.02-10-4

qo 0.120 0.365 0.250 0.201

E 252 217 1.910-1 3.110-6

Для тестирования работы алгоритма была решена прямая задача (22) для набора параметров а=220, Ь=60, Б=1.0-10-4, Чо=0.20 и 1=1, Т=1. Результат решения для расхода и концентрации мобильной компоненты в точке пространства х=Ь был сохранен с шагом по времени А/=1/850. Таким образом, полученное решение представлено в виде набора из 850 точек. Полученные данные были использованы в качестве эталонного набора данных (заменяющих предполагаемые экспериментальные данные). Эти данные графически представлены на рис. 5, 6 сплошной линией.

Тестирование алгоритма проводилось для начального набора параметров а=212, Ь=80, Б=6.710-4, Чо=0.12. Оно показало, что после 40 итераций алгоритм сошелся к искомому набору параметров, характеризующих набор эталонных данных. Эволюция значения целевой функции во время работы алгоритма представлена на рис. 3.

Из рис. 3 видно, что при заданных начальных параметрах значение целевой функции Е ~ 103, после трех итераций достигается локальный минимум, в котором значение целевой функции становится порядка Е ~ 1. После п=13 итерации функция имеет скачок. Скачок значения целевой функции вызван большим шагом в пространстве параметров, необходимым для выхода из окрестности локального минимума. Начиная с п=17 итерации значение целевой функции монотонно уменьшается до итерации п=24, где Е ~ 10-5. Исходя из последующей работы алгоритма можно предположить, что это - глобальный минимум, так как целевая функция не изменяется вплоть до завершения работы алгоритма при п=40 итерации. Изменение значений параметров в ходе работы алгоритма представлено в таблице. Можно видеть, что при итерации п=40 значения параметров фактически совпадают со значениями, выбранными для эталонных данных. Максимальная разница между ними составляет Л-7-10-2, что свидетельствует о сходимости и эффективности алгоритма.

8. Заключение

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

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

Исследование выполнено за счет гранта Российского научного фонда (проект № 20-11-20125).

Список литературы

1. Колчанов Н. В., Евграфова А. В., Казакова Е. А., Колчанова Е. А. Влияние примесей NaCl на расход протекающей сквозь пористую среду воды // XXII Зимняя школа по механике сплошных сред Пермь, 22-26 марта 2021 г.: Тез. докл. (Электронный ресурс). Пермь, ПФИЦ УрО РАН, 2021, С. 170. (дата обращения: 26.03.2021)

2. Bromly M., Hinz C. Non-Fickian transport in homogeneous unsaturated repacked sand // Water Resources Research. 2004. Vol. 40. No. 7. WR002579. DOI: 10.1029/2003WR002579

3. Chavent G., Dupuy M., Lemmonier P. History matching by use of optimal theory // Society of Petroleum Engineers Journal. 1975. Vol. 15. N. 01. P. 74-86. DOI: 10.2118/4627-PA

4. Гилл Ф., Мюррей У., Райт М. Практическая оптимизация. M.: Мир, 1985. 509 c.

5. Акулич И. Л. Математическое программирование в примерах и задачах. M.: Лань, 2009. 347 c.

6. Fletcher R., Reeves C. M. Function minimization by conjugate gradients // The Computer Journal. 1964. Vol. 7. N. 2. P. 149-154. DOI: 10.1093/comjnl/7.2.149

7. Fletcher R. Practical methods of optimization. Wiley, 2013. 456 p.

8. Nocedal J., Wright S. Numerical optimization. Springer Science & Business Media, 2006. 664 c.

9. Chavent G., Lemonnier P. Identification de la non-linéarité d'une équation parabolique quasilinéaire // Applied mathematics and Optimization. 1974. Vol. 1. N. 2. P. 121-162. DOI: 10.1007/BF01449027

10. Будак Б. М, Гапоненко Ю. Л., Малышева Г. Ю., Рубан П. И. Об одном методе решения экстремальных задач с ограничениями на фазовые координаты // Журнал вычислительной математики и математической физики. 1974. Т. 14. Вып. 3. С. 779-783.

11. Deans H. A. A mathematical model for dispersion in the direction of flow in porous media // Society

of Petroleum Engineers Journal. 1963. Vol. 3. N. 01. P. 49-52. DOI : 10.2118/493-PA

12. Van Genuchten M. T., Wierenga P. J. Mass transfer studies in sorbing porous media I. Analytical solutions 1 // Soil science society of america journal. 1976. Vol. 40. N. 4. P. 473-480. DOI: 10.2136/sssaj1976.03615995004000040011x

13. Schumer R., Benson D. A., Meerschaert M. M., Baeumer B. Fractal mobile/immobile solute transport // Water Resources Research. 2003. Vol. 39. N. 10. DOI : 10.1029/2003WR002141

14. Gouze P., Le Borgne T., Leprovost R., Lods G., Poidras T., Pezard P. Non-Fickian dispersion in porous media: 1. Multiscale measurements using single-well injection withdrawal tracer tests //Water Resources Research. 2008. Vol. 44. N. 6. DOI: 10.1029/2007WR006278

15. Selim H. M. Prediction of contaminant retention and transport in soils using kinetic multireaction models // Environmental health perspectives. 1989. Vol. 83. P. 69-75. DOI: 10.2136/sssaj2006.0422

16. Maryshev B. S. The linear stability of vertical mixture seepage into the close porous filter with clogging // Fluid Dynamics Research. 2016. Vol. 49. N. 1. P. 015501. DOI: 10.1088/01695983/49/1/015501

17. Roth E.J., Gilbert B., Mays D.C. Colloid deposit morphology and clogging in porous media: fundamental insights through investigation of deposit fractal dimension // Environmental Science & Technology. 2015. Vol. 49. P. 12263. DOI: 10.1063/1.4935576

18. Carman P. C. Fluid flow through granular beds //Trans. Inst. Chem. Eng. 1937. Vol. 15. P. 150166.

19. Kozeny J. Uber kapillare leitung der wasser in boden // Royal Academy of Science, Vienna, Proc. Class I. 1927. Vol. 136. P. 271-306.

20. Darcy H. Les fontaines publiques de la ville de Dijon: exposition et application.. Victor Dalmont, 1856.

21. Langmuir I. The adsorption of gases on plane surfaces of glass, mica and platinum // Journal of the American Chemical society. 1918. Vol. 40. N. 9. P. 1361-1403.

22. Самарский А. А. Введение в численные методы. М.:.Лань, 2009.

References

1. Kolchanov N.V., Evgrafova A.V., Kazakova E.A., Kolchanova E.A. Vliyanie primesey NaCl na ra-shod protekayushey skvoz poristuyu sredu vody (The effect of NaCl admixture to filtration flux of water through porous medium) Proc. All-Russian Conf. "Zimnyaya shkola po mehanike sploshnyh sred" (Winter school on continuous media mechanics) 2021, Perm, Russia, Perm Federal Re-

search Center of the Ural Branch of the Russian Academy of Sciences, p. 170.

2. Bromly M., Hinz C. Non-Fickian transport in homogeneous unsaturated repacked sand. Water Resources Research, 2004, vol. 40, no. 7, WR002579, DOI: 10.1029/2003WR002579

3. Chavent G., Dupuy M., Lemmonier P. History matching by use of optimal theory Society of Petroleum Engineers Journal, 1975, vol. 15, no. 1, pp. 74-86. DOI: 10.2118/4627-PA

4. Gill P. E., Murray W., Wright M. H. Practical optimization. Philadelphia, USA, SIAM pub., 2019, 401 p.

5. Akulich I. L. Matematicheskoe programirovanie v primerah i zadachah (Mathematical programming: problems and exercises), Moscow, Russia, Lan', 2009, 347 p.

6. Fletcher R., Reeves C. M. Function minimization by conjugate gradients The computer journal., 1964, vol. 7, no. 2, pp. 149-154. DOI: 10.1093/comjnl/7.2.149

7. Fletcher R. Practical methods of optimization. , John Wiley & Sons, 2013, 456 p.

8. Nocedal J., Wright S. Numerical optimization,. Springer Science & Business Media, 2006, 664 p.

9. Chavent G., Lemonnier P. Identification de la non-linéarité d'une équation parabolique quasilinéaire Applied mathematics and Optimization, 1974, vol. 1, no. 2, pp. 121-162. DOI: 10.1007/BF01449027

10. Budak, B. M., Gaponenko, Y. L., Malysheva, G. Y., & Ruban, P. I. A method of solving extremal problems with constraints on the phase coordinates. USSR Computational Mathematics and Mathematical Physics, 1974, vol. 14, no 3, pp. 239-243.

11. Deans H. A. A mathematical model for dispersion in the direction of flow in porous media Society of Petroleum Engineers Journal, 1963, vol. 3, no. 1. pp. 49-52. DOI : 10.2118/493-PA

12. Van Genuchten M. T., Wierenga P. J. Mass transfer studies in sorbing porous media I. Analytical solutions 1 Soil science society of America journal, 1976, vol. 40, no. 4, pp. 473-480. DOI: 10.2136/sssaj1976.03615995004000040011x

13. Schumer R., Benson D. A., Meerschaert M. M., Baeumer B. Fractal mobile/immobile solute transport Water Resources Research, 2003, vol. 39, no. 10, WR002141. DOI : 10.1029/2003WR002141

14. Gouze P., Le Borgne T., Leprovost R., Lods G., Poidras T., Pezard P. Non-Fickian dispersion in porous media: 1. Multiscale measurements using single-well injection withdrawal tracer tests Water Resources Research. 2008, vol. 44, no. 6. DOI: 10.1029/2007WR006278

15. Selim H. M. Prediction of contaminant retention and transport in soils using kinetic multireaction models Environmental health perspectives, 1989, vol. 83, pp. 69-75. DOI: 10.2136/sssaj2006.0422

16. Maryshev B. S. The linear stability of vertical mixture seepage into the close porous filter with clogging Fluid Dynamics Research. 2016, vol. 49, no. 1, 015501. DOI: 10.1088/01695983/49/1/015501

17. Roth E.J., Gilbert B., Mays D.C. Colloid deposit morphology and clogging in porous media: fundamental insights through investigation of deposit fractal dimension Environmental Science & Technology. 2015, vol. 49, pp. 12263. DOI: 10.1063/1.4935576

18. Carman P. C. Fluid flow through granular beds Trans. Inst. Chem. Eng. 1937, vol. 15, pp. 150166.

19. Kozeny J. Uber kapillare leitung der wasser in boden Royal Academy of Science, Vienna, Proc. Class I, 1927, vol. 136, pp. 271-306.

20. Darcy H. Les fontaines publiques de la ville de Dijon: exposition et application. Victor Dalmont, Paris, 1856, 649 p.

21. Langmuir I. The adsorption of gases on plane surfaces of glass, mica and platinum Journal of the American Chemical society, 1918, vol. 40, no. 9. pp. 1361-1403.

22. Samarsky A. A. Vvedenye v chislennye metody (Intoduction to numerical methods), Moskow, Lan', 2009, 288 p. (In Russian).

Просьба ссылаться на эту статью в русскоязычных источниках следующим образом: Хабин М. Р., Марышев Б. С. Идентификация параметров транспорта примеси через пористую среду с учётом закупорки // Вестник Пермского университета. Физика. 2021. № 3. С. 44-55. doi: 10.17072/19943598-2021-3-44-55

Please cite this article in English as:

Khabin M.R., Maryshev B.S.. Identification of the parameters of transport through porous media with clogging . Bulletin of Perm University. Physics, 2021, no. 3, pp. 44-55. doi: 10.17072/1994-3598-2021-3-44-55

Сведения об авторах

1. Михаил Романович Хабин, студент, инженер кафедры теоретической физики, Пермский государственный национальный исследовательский университет, ул. Букирева, 15, Пермь, 614990

2. Борис Сергеевич Марышев, канд. физ-мат. наук, н.с., Институт механики сплошных сред УрО РАН, Пермь, ул. Академика Королева, 1, Пермь, 614013; доцент кафедры теоретической физики, Пермский государственный национальный исследовательский университет, ул. Букирева, 15, Пермь, 614990

Author information

1. Mikhail R. Khabin, student, engineer at the Theoretical Physics Department, Perm State University, Bukirev str. 15, Perm, Russia, 614990

2. Boris S. Marychev, Candidate of Physical and Mathematical Sciences, Researcher, Institute of Continuous Media Mechanics UB RAS, Akademika Koroleva str. 1, Perm, Russia, 614013; Professor Associate at the Theoretical Physics Department, Perm State University, Bukirev str. 15, Perm, Russia, 614990