О методах идентификации водопроводимости пластов Текст научной статьи по специальности «Математика»

Математика. Физика

УДК 532. 546+518.5

О МЕТОДАХ ИДЕНТИФИКАЦИИ ВОДОПРОВОДИМОСТИ ПЛАСТОВ М.У. Мурзакматов

Кафедра «Прикладная математика»

Иссык-Кульский государственный университет им. К. Тыныстанова, Республика Кыргызстан

Представлена членом редколлегии профессором В.И. Коноваловым

Ключевые слова и фразы: водоносный пласт; водопроводимость; гидрогеологические параметры; подземные воды; пористая среда; фильтрация.

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

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

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

В зависимости от степени изученности области фильтрации в некоторых точках могут задаваться значения напорной функции Н(х, у), коэффициента фильтрации к(х, у), водопроводимости Т(х, у), упругой водоотдачи |т(х, у) и инфильтрации /(х, у). По этим неполным данным необходимо восстановить картину течения жидкости, то есть идентифицировать фильтрационный поток.

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

Н(х1, у) = Н, (х, у )е Кд, 1 = l, 2, ••• , Ч, (1)

Т(х1, у) = Т'Э, (х1, у )е ^, 1 = 1 2, - , Р , (2)

ггЭ Т^Э __

где Н1 , Т' - значения указанных величин, полученные из эксперимента или путем наблюдений; Кд и Sp - некоторые дискретные множества, состоящие из д и

р точек, соответственно. Требуется доопределить значения напоров и водопрово-димости в области фильтрации Б, где заданы множества Кд и Sp , считая остальные параметры потока известными. Такая задача рассматривалась в работе [1].

В стационаром случае плановый фильтрационный поток описывается уравнением

ЬН = /, (х, у) е Б (3)

с краевым условием

Н = а, (х, у) е Г = дБ. (4)

Здесь

ь = _А{ТА] _А(ТА]+е, I = Тд_р,

дх ^ дх) ду { ду) ^ дп н

Н = Н(х, у) - функция напора, м; Т = Т(х, у) - водопроводимость водоносного пласта, м2/сут; <2 = 2(х,у) - функция, учитывающая переток из ниже- и вышележащих пластов, м/сут;/=/ (х, у) - функция источников и стоков, м/сут; а = а(х, у), м2/сут

и Р = Р(х, у), м/сут - заданные функции; Б - область фильтрации в плане; Г - ее

граница; д/дп - производная по внешней нормали к границе области.

Задачу (3), (4) при известных функциях Т, 2, / Р, а и известной границе Г можно решить методом конечных элементов [4]. Плоскую область Б произвольным образом разбиваем на т треугольных элементов. Пусть элемент Ае имеет

своими вершинами точки 1, у, к с координатами (х', у'), (ху, у у), (хк, ук).

Внутри этого элемента решение задачи находится в виде линейной комбинации

Не (х, у) = N (х, у)Н + (х, у)Ну + Ык (х, у)Нк , (5)

где Н8 = Н(х^, у!.), 5 = 1, у, к ,

N = а1 + Ъ1х + с1у; N у = а у + Ъух + Суу; Ык = йк + Ькх + Ску;

хуук _ хкуу , у} _ ук хк _ ху

а1 = —-----------; Ъ1 = —-------; с =-----------; ...;

1 2Ае 1 2Ае 1 2Ае

остальные коэффициенты выписываются с помощью круговой постановки индексов 1, у, к; Ае - площадь треугольника Ае .

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

п

Н (х, у) » Нп (х, у) = X НК1 (х, у) ® . (6)

1=1

Подставим в задаче (3), (4) вместо Н функцию Нп (х, у) из (6), применяем к этой задаче обобщенный принцип Галеркина и после упрощения с помощью формулы Грина получаем систему линейных алгебраических уравнений относительно

Н,:

п

Еа]1Н1 =ъ]; , = l,2,...,n, (7)

1=1

где

«„ = JJ Ш N, ^ )ds + JJNjNjQds + J NjNjbds ;

D D Г

bj = JJ fNjds + J aNjds;

, > +

В Г

ЭМ ЭМ,- ЭМ ЭМ,-д( М , М,- ) = + ЭМ^_± .

^ Эх Эх Эу Эу

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

К задаче идентификации водопроводимости можно использовать различные подходы. Мы рассмотрим три из них.

1° Этот подход условно назовем методом итераций, хотя итерации используются и в других методах. Алгоритм этого метода предложен в работе [1], который состоит в последовательном решении задачи (3), (4) относительно Н(х, у) и Т(х, у) с использованием условий (1) и (2). Алгоритм решения задачи строится следующим образом.

В сеточной области функцию Т(х, у) представим в виде

п

Т(х, у) » Тп(х, у) = Ет,М, (х, у) (8)

1=1

и в первом приближении задаем ее значения произвольным образом из интервала

а в узлах, совпадающих с точками множества , задаем значения

Тэ . Применяя к задаче (3), (4) с условием (1) обобщенный принцип Галеркина,

получаем систему уравнений

я

JJМ1(ЬН-/№+Е Н -нэ + JМ1(1Н-а)й& = 0, I = 1, 2, ..., п. (9)

в ,=1 Г

После подстановки вместо Н функции Нп (х, у) из (6) с учетом формулы Грина приходим к системе

п

Еа1]Н] = 81, 1=1 2 .., п (10)

,=1

где

atJ = JJ Tq( Nt, Nj )ds + JJNtNjQds- J NtNjbds + dj;

D D Г

qt

gt = JJNtfds+ JNtads + ^dsH3s .

D

s=1

Здесь - число узлов шаблона, центром которого является узел 1 (включая сам узел 1);

§s =

1, если значение Н, задано,

0, если Нэ не задано, 5 = 1, 2, ..., п.

Решение Н(1)(х, у) системы (10) является первым приближением напорной

функции Н(х, у). Используя полученные значения Н(1 (1 = 1, 2, ..., п) и условие

(2), решаем задачу (3), (4) относительно водопроводимости Т(х, у). Здесь, кроме выполнения внутренних условий (2), потребуем также, чтобы вариация искомой функции стремилась к нулю. Это условие должно обеспечить требуемую гладкость функции Т(х, у). Тогда вместо (9) получаем систему

Р

JJ М1 (ЬНп - /)ё<5 + JJ М1 ЪТёО + Е Т5 - тэ + J М1 (1Нп - а)Ж = 0, 1 = 1, 2, ..., п,

D

D

s=1

Г

где ЪТ = Т - Т, Т - значения функции Т, полученные из предыдущей итерации. Применяя формулу Грина и подставляя вместо Тп ее разложение (8), приходим к системе

п

Е с1]Т]=, 1=1 ^ ..., n,

,=1

(11)

где

ctJ = JJ NJq( Nt, Hn )d s + JJ NtNJd s + gJ;

DD

dt = JJNfds- JJNtHnQds + JJNtTnds + gT + J Nt ($H„ + a)ds;

(12)

0 0 0 Г

11, если значение Т/ задано,

15 = 1

10, если значение Т, не задано, 5 = 1, 2, ..., р.

Решив систему (11), получаем второе приближение Т(2)(х, у) и подставляя его вместо Т, снова приходим к системе (10). Решение этой системы, т.е. второе приближение функции Н(х, у) используем для нахождения Т(3) (х, у) из системы (11) и т.д. Итерационный процесс продолжается до выполнения условия

тах

< е, 1 = 1, 2, ..., п,

где V - номер итерации; е > 0 - заданное малое число.

(13)

Г

Наибольшую трудность в этой задаче представляет решение системы (11). Как видно из формул (12), коэффициенты этой системы мало отличаются друг от друга, вследствие чего ее матрица является почти вырожденной. Для решения такой плохо обусловленной системы необходимо применять метод сингулярного разложения матрицы [6].

2° Этот подход основан на методе регуляризации А.Н. Тихонова [5]. Задача сводится к нахождению функции Т(х, у) из уравнений (3), (4), сообщающей в области Б минимум функционалу [3]

q г

Ф(Т) = £[ 2=1

Hi (T) - H

Е(

j=1

Tj Tj

)2 +71ST |2.

(14)

Здесь НI (Т) - расчетные значения напоров, которые находятся как решение задачи (3), (4); у - параметр регуляризации.

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

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

~ П ~ 7Г

H(T) = H(T) + Е (T - Ts) — + R2 (AT);

s=1 dTs

(15)

где Т8 -значение функции Т в точке (хда у), полученное в предыдущей итерации; ^2(ДТ) - остаточный член разложения. Подставляя (15) в (14) и используя необходимое условие минимума функции многих переменных, получаем систему линейных алгебраических уравнений

ЭФ(Т)

dTk

'l

=Е

i=1

~ П ~ 7Г

нi + Е (Ts - Ts)-------------------------l- - Hi

2 ^ s dTs 1

k = 1, 2,

s=1

ъH ~

l +mk (Tk - Tk)+g(Tk - Tk) = o,

dTk

с коэффициентами

ЕcksTs = dk , k = 1 2, -,

(16)

s=1

i=1 V *Tk

+g+mk

и с правыми частями

i=1

dHi

i dTs

s=1 s J

-H~+mkTk+g Tk ■ °Tk

Здесь

mk=<

11, если Tэ задано,

I 0, если Tэ не задано.

2

n

n

Г

q

n

Вычислительная процедура осуществляется в следующем порядке. Используя начальные значения Н^п £ Н;-0) £ НШах, Тпт £ Т'^ £ ТШах ( = 1, 2, ..., п) в качестве нулевого приближения, решается задача (3), (4) и определяется первое приближение Н(1). Затем, придавая приращение АТ функции Т(х, у), находим второе приближение Н(2). Это дает возможность приближенно определить производную дН/дТ и решить систему (16) при некотором значении параметра регуляризации у. Найденные значения Т(х, у) используются для определения следующего приближения Н. Этот процесс продолжается до установления фильтрационного потока, т.е. до выполнения условия (13). Если при этом полученные значения напоров в пределах допустимых погрешностей не совпадут с данными

экспериментальными значениями Нэ, то итерация продолжается по параметру регуляризации у, выбор которого может быть осуществлен методом невязок [5].

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

Согласно [2], идентификация водопроводимости Т(х, у) проводится в следующей последовательности.

Образуем начальное приближение искомой функции Тг(0), где

ттТэ £ Т£ тахТгэ, г = 1, 2, ..., п и составим функцию ТП°\х, у) по форму-

г г

ле (8). Решив краевую задачу (3), (4) с функцией ТП°\х, у), находим начальное

приближение напорной функции НП0) (х, у).

Пусть теперь правые части уравнения (3) и краевого условия (4) получают малые возмущения 8/ и 8а, соответственно. Тогда задача (3), (4) переходит к задаче

ЬН' = /'; Н' = Н + 8Н; /' = / + 8/; (х, у) е Б; (17)

1Н' = а ; а' = а + 8а; (х, у) е Г. (18)

Из (17) и (18) с учетом (3) и (4) приходим к краевой задаче относительно вариации 8Н :

Ь8Н = 8/; (х, у) е Б; 18Н = 8а; (х, у) е Г. (19)

Задача (19) при известных вариациях 8/ и 8а позволяет построить поле вариаций напоров 8Н.

Решив задачу (19) при той же функции ТП°^(х, у) и заданных функциях

8/(х, у) и 8а(х, у), находим вариацию напорной функции 8Н(х, у) и преобразуем в функцию

Н'(х, у) = Н(0) (х, у) + 8Н (х, у).

Теперь рассмотрим случай, когда малые возмущения получают не только правые части уравнений (3) и (4), но и их коэффициенты Т(х, у), <2(х, у) и Р(х, у).

Это равносильно изменению гидрогеологических характеристик пласта. Пусть вместо задачи (3), (4) имеем краевую задачу

где

Ь И = /"; Ь = Ь+ 8Ь; Г = /' + 8/'; (х, у) є Б; 1'И' = а ; Ґ = I + 81; а=а+8а ; (х, у) є Г,

8Ь = -—\ 8Т—1 \ 5Т — 1 + 80, 81 = 8Т—-8В.

дх ^ дх) ду \ ду) дп н

Из (20), (21) с учетом (17), (18) получаем краевую задачу 8ЬИ' = 8/'; (х, у) є Б; 8И' = 8а'; (х, у) є Г.

(20)

(21)

(22)

Поскольку функция Н (х, у) известна, отсюда получаем уравнение относительно вариации искомой функции 8Т(х, у) . Решение задачи (22) ищем в виде

п

8Т(х, у) » 8Тп(х, у) = £ 8Т]Ы] (х, у). (23)

] =1

Применяя обобщенный принцип Г алеркина к задаче (22), имеем

ЛЫ,(8ЬН'-8/)йо + |Ы,(81Н'-8ай = 0, г = 1, 2, ..., п.

Б Г

Подставляя вместо 8Т ее разложение (23) и используя формулу Грина, приходим к системе линейных алгебраических уравнений относительно 8Тг

£ Л Ыуд( Ы,, Н' )й о 8Т} = Л Ы, 8/ й о - Ц Ы,Н' 8<2 й о-1Ы, (Н '8р - 8а )&, (24) У=1 V Б У Б Б Г

I =1, 2, ..., п .

Система (24) является плохо обусловленной, так как функция д(Ы,, Н') не

обеспечивает системе диагонального преобладания и ее определитель близок к

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

п ( 1 Р

£ ЦЫ]д(Ы1,Н')йо 8Т] +т||Ы,8Тпйо + у£ (8Т, - Т,3 + Т,) =

У=1 V Б У Б

s=1

Л N. df'ds - JJ N. H'dQd s - J N. (Hdp - da')ds,

i = 1,2,..., n,

D

D

или

где

j aij dTj = fi , .= 1, 2,■■■, n, j=1

aij = JJ Njq( Nj, H' )d s + mJJ N.Njd s + gj

(25)

(2б)

D

D

Ы,8/йо-ЦЫ,Нъдйо-1Ы, (Н'8Р-8а')йя + у, (Т,э -Т,),

Б Б Г

П, если Т, задано,

У, = 1 э 5 = >, J, К

[0, если Т, не задано,

т - положительное число, выполняющее роль параметра регуляризации; Т - значение Т из предыдущей итерации.

Работа алгоритмов всех трех методов апробирована на серии тестовых при-

2 2

меров. Рассмотрим один из них. В круге х + у £ 0, 25 решается задача (3), (4) с внутренними условиями (1), (2) (в методе возмущений условие (1) вообще не используется). Круг разделен на 54 элемента, число узлов сетки - 37, из них 18 -граничных. Элементами являются равносторонние треугольники со стороной,

2 2

равной 0,2. В этой области заданы функции Н(х, у) = х + у + 5, 2(х, у) = 0,

/ (х, у) = -40(2х + 2у +1). Искомой функцией является Т(х, у) = 10( х + у +1). Числа д и р в условиях (1) и (2) (т.е. число узлов, в которых известны экспериментальные (точные) значения функций Н и Т соответственно) задавались в различных сочетаниях (табл. 1).

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

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

Таблица 1

Относительные погрешности в определении водопроводимости различными методами (в процентах)

q P Метод итераций Метод регуляризации Метод возмущений

1? 1? 8,8 б,3 ?,1

1? 13 1О,2 9,б ?,4

1З 1? 11,О 1О,4 ?,З

1З ? 11,? 11,2 8,5

? 13 9,? 1О,б ?,?

? 4 1?,4 12,3 8,б

4 ? 11,5 9,1 8,?

4 4 23,5 12,? 9,3

= JJ

1 Джаныбеков, Ч. Моделирование гидрогеодинамических процессов с применением ЭВМ / Ч. Джаныбеков. - Фрунзе : Илим, 1989. - 183 с.

2 Джаныбеков, Ч. Методы идентификации гидрогеологических параметров и прогнозирования процессов загрязнения подземных вод / Ч. Джаныбеков, М.У. Мурзакматов. - Бишкек : Илим, 2005. - 180 с.

3 Мурзакматов, М. У. Идентификация водопроводимости пласта методом регуляризации / М.У. Мурзакматов, Ж.М. Мамыров // Докл. НАН Республики Казахстан. - 2004. - №6. - С. 82-87.

4 Сегерлинд, Л. Применение метода конечных элементов / Л. Сегерлинд. -М. : Мир, 1979. - 392 с.

5 Тихонов, А.Н. Методы решения некорректных задач / А.Н. Тихонов, В.Я. Арсенин. - М. : Наука, 1974. - 233 с.

6 Форсайт, Дж. Машинные методы математических вычислений / Дж. Форсайт, М. Малькольм, К. Моулер. - М. : Мир, 1980. - 279 с.

To the Question of Identification Methods of Aquifers Permeability

M.U. Murzakmatov

Issykkul State University after K. Tynynstanov, Kyrgyz Republic

Key words and phrases: aquifer; filtration; ground water; hydro-geological parameters; permeability; porous media.

Abstract: The algorithms of three methods of identification of aquifers permeability are given. Results of solution to the test task are compared.

Uber die Methoden der Identifizierung der Wasserleitfahigkeit der Schichten

Zusammenfassung: Es werden die Algorithmen der drei Methoden der Identifizierung der Wasserleitfahigkeir der Wasserhorizonte angefuhrt und es werden die Ergebnisse des Beschlusses der Prufungsaufgabe verglichen.

Sur les methodes de l’identification de la permeabilite des aquiferes

Resume: Sont cites les algorithmes des trois methodes de l’identification des nappes aquiferes et sont compares les resultats de la solution du probleme de test.