CC BY
16НАУКА И ОБРАЗОВАНИЕ, 2015, №1
УДК 532.5 (075)
Методы идентификации нестационарных математических моделей
теории фильтрации
Э.А. Бондарев, А.Ф. Воеводин*
Институт проблем нефти и газа СО РАН, г. Якутск *Институт гидродинамики СО РАН, г. Новосибирск
Математическими моделями движения жидкостей и конвективной диффузии примесей в пористых средах являются дифференциальные уравнения эволюционного типа. Для определения их параметров по дополнительной информации о поведении искомых функций (данные измерений) предложено два метода, один из которых применим для определения постоянных параметров, а второй -для параметров, зависящих от пространственных координат. На конкретных примерах показано, что эффективность методов существенно возрастает с увеличением числа пространственных измерений, тогда как увеличение их частоты по времени ведет к повышению чувствительности - процесс восстановления приобретает колебательный характер.
Ключевые слова: математическое моделирование, методы идентификации, теория фильтрации.
Mathematical models offluid flow and contaminations spreading in porous media are partial differential equations of evolution type. Two methods of their parameters identification have been proposed: the first one for constant parameter values, the second - for distributed ones. In terms of boundary value problems for the equation of transient fluid flow in porous media it has been shown that efficiency of the methods strongly depends on the number of measuring transducers whereas sample rate influences the methods' response.
Key words: mathematical models, identification methods, filtration theory.
Гидравлические модели, соответствующие нестационарным режимам течения, описываются дифференциальными уравнениями в частных производных эволюционного типа. Задачи идентификации параметров таких моделей начали рассматриваться сравнительно недавно в силу значительных вычислительных трудностей и отсутствия общей теории для таких задач. Используемые в настоящее время методы их решения и в идейном плане, и в техническом осуществлении сходны с методами идентификации моделей, описываемых обыкновенными дифференциальными уравнениями [1].
Пусть общая эволюционная задача имеет вид ду
dt
= f fc - У, У^ А)
(1)
с необходимыми начальными и краевыми условиями в области t > 0, 0 < x < 1. Здесь y - вектор искомых функций, t - время, x - пространственная переменная; нижние индексы обозначают частные производные первого и второго порядка.
Требуется определить вектор неизвестных параметров А = const из условия минимума функционала
БОНДАРЕВ Эдуард Антонович - д.т.н., проф., г.н.с., [email protected]; *ВОЕВОДИН Анатолий Федорович - д.ф.-м.н., проф., г.н.с., [email protected] .
1 1 t J=JJi y - -У) G(y ~ y)dxdt-
(2)
При известном начальном приближении Л" и соответствующем ему у" исходную задачу (1) можно свести к последовательности линейных задач
Ь/+1 = а^,х)Л+1 + Ь(1х), (3) где L - линейный дифференциальный оператор, имеющий вид
L =
д_ dt'
fy
dx2
fy
Ух
д_
dx
fs
J y
а ait, x) и bit, x) - соответственно матричная и
вектоРная функции a = f, b = f s - fSy y Xx -
S Г S Г OS A
ухУх - Jyy - JaÀ . Аналогично линеаризуются граничные условия.
Решение (3), используя идеи метода матричной прогонки по параметру [1], представим в виде
у+1 = х + ^ (4)
В результате из (3) получим для определения X и z две линейные краевые задачи:
LX = aix) ; (5)
Lz = bix). (6)
Решение задач (5), (6) и формула (4) позволяют в явном виде получить зависимость
s
J (х+1) и в дальнейшем определить Х+1 из необходимого условия стационарности по формуле
X's+1 =
О 0
T 1 1 T 1 ,
Л X'GXdxdt JJ(y - z) GXdxdt.
О О
На этом заканчивается один шаг итерационного процесса. Следующие приближения Х+2, Л"+3 и т.д. определяются аналогично.
Этот вариант метода идентификации будет нами использоваться при решении некоторых конкретных задач. Основное его достоинство заключается в том, что на каждом итерационном шаге решаются краевые задачи (5) и (6), отличающиеся только правой частью, при этом если исходная задача при заданных X корректно разрешима, то и задачи (5), (6) корректно разрешимы.
Для восстановления распределенных, т.е. зависящих от пространственной переменной, параметров модели в работах [2, 3] предложен градиентный метод, который базируется на необходимых условиях оптимальности (равенство нулю первой вариации функционала). Пусть модель некоторого процесса описывается уравнением
дУ = /(^х,у,ух,у„, Л(х}лх (х)), (7)
дt
где неизвестный параметр есть функция
переменной х . Начальные и граничные условия заданы в виде:
у(х,0) = у0 (х); (8)
§0 у,ух ) = 0 при X = 0; (9)
gl (г,у,ух )= 0 при х = 1. (10)
Требуется выбрать
Л(х) таким образом, чтобы минимизировать функционал (2). Применяя стандартную технику [3], можно показать, что первая вариация функционала имеет вид
j=J iJ[v ¡к)
+
fx5(x -1)-( fx5(x )]йф dx,
где ¿(х) - дельта-функция Дирака; у - сопряженная вектор-функция, удовлетворяющая уравнению
дУ = -у'л Ау'/„)х )„ (11)
и условиям
v(x,T ) = 0;
( fyx fyx )x fygyxg0y = 0
при x = 0;
(12)
(13)
V
' fy-(v fy„ 1 -V fy xS-lSly = 0
при х = 1. (14)
Таким образом, для того чтобы приращение ¿И было отрицательно на каждой итерации, достаточно задать приращение вектора параметров в виде
т
5Л(х ) = -W (x)j(v /я-((' К )х
0
+ V' fK8(x - 1)-v' fx8(x )}lt\,
(15)
где W (х) - положительно-определенная матрица весомых множителей размерности г X г .
Численный алгоритм нахождения Х(х ) можно построить следующим образом:
1) выбираем начальное приближение Л° (х) и матрицу весовых множителей W (х);
2) решаем систему (7)-(10) от t = 0 до £ = Т . Находим значение функционала J(X (х)). Решаем систему (11)—(14) от t = Т до t = 0;
3) по формуле (15) вычисляем 8Х (х) и (х) = Х (х) + 5А? (х);
4) повторяем шаг 2 до выполнения условия
- Г )/ Г <£.
Дальнейшее изложение ведется на конкретном примере модельного уравнения теории фильтрации с постоянными коэффициентами. Это уравнение соответствует линеаризованному уравнению Буссинеска для прогноза уровня грунтовых вод или уравнению упругого режима фильтрации.
В области х &\0,1 ], t &\0,Т ] рассмотрим уравнение
дп „ д2u r! \
--К—- = f (x,t),
dt dx2 V У
(16)
в котором коэффициент X = const подлежит определению.
Пусть x1,x2,...,xn - характерные точки мелиорируемой территории, в которых производится измерение уровня грунтовых вод. Данные этих измерений - известные функции времени Z, (t) 1 < j < n.
Как и прежде, будем искать параметр К из условия минимума функционала
n т
J(A=1LJ¥Xj,S)-(¿)!<%. (17)
0
Кроме того, для уравнения (16) обычным образом задаем начальные и граничные условия.
Условие стационарности функционала (17) будет иметь вид
сСЛ (Л)_ СЛ
п Т
¡Их,,?)—г, (4)]
(X
, =1 о '
м{х], = о.
(18)
Здесь, как и ранее, мы ввели специальное обозначение для производной функции состояния и по параметру Л
V = иЛ. (19)
Разложим в ряд функцию и в окрестности Л с точностью до членов второго порядка
из+1 ) = и (х,г)+(Л+1 — Л (х,г). (20)
Для сокращения записи здесь и далее считается, что нижний индекс 5 у функции означает, что она вычисляется при значении Л = Л .
Подставляя разложение (20) в соотношение (18), линеаризуем это соотношение относительно параметра Л1 следующим образом:
2Е/к 4)+(Л+1 -Л к X, 4)
1 0 — z,.
(ф* (х;, 4^4 = 0.
Отсюда легко вычислить следующее приближение Л+1, если функции и,, (л, I) и ws (к, t) известны
X (х],4)—щ (хр4)+
Л+1 =
] 0
+ (4%>* (х,,4^4
Т ■
XI ^ (х],4)44
(21)
] о
Чтобы получить уравнение для функции ws (х, t), проведем линеаризацию исходного уравнения (16) относительно решения на нижнем итерационном слое:
ди,+1 д2и* д2и* *+1 ■ = Л —^ + —* 8Л +
ы
+ Л
*+1
дх2 дх2
( д 2и
*+1
дх2
д 2п, ~дх2
\
(22)
+1,
где ¿Л = Л —Л.
Подставляя в уравнение (22) разложение (20) и приравнивая коэффициенты при ¿Л нулю, получим следующую систему уравнений:
^ = Л д+
дг дг
= Л
дх2 д 2~м„ д2и„
дх2
■ + -
дх2
(23)
(24)
Уравнения (23) и (24) при заданном значении Л могут быть решены последовательно одним из известных численных методов. Граничные и начальные условия для функции м могут быть получены из соответствующих условий для функции и путем дифференцирования их по параметру Л.
Численную реализацию изложенного метода рассмотрим на примере определения параметра Л в однородном уравнении (16)
ди д ^и
-= Л-0 < X < 1
д дх2
с начальным условием и(х, 0) = 0 и граничными
условиями ди = о при х = 0, и = 1 при х = 1.
дх
Эта задача имеет аналитическое решение
( \ , 4 1) и(х, t) = 1--X -—— ехр
п i=1 21 +1
Л(21 + l)2п2t
х cos
(21 +1). 2
пх
которое использовалось для вычисления «данных измерений» г(Х], ¿к) при Л = 0,2 . Эти величины вычислялись в пяти точках пространства х. = 0; 0,2; 0,4; 0,6; 0,8 и в дискретные моменты времени tk , число которых варьировалось. Кроме этого варьировалась величина временного интервала Т, на котором осуществлялась идентификация. Как показали вычисления, даже при задании начального приближения Л0 в точности равного 0,2, рассмотренный алгоритм расходится. По мере увеличения количества итераций параметр Л колеблется возле точного значения с возрастающей амплитудой. Поэтому для решения данной задачи нами был использован модифицированный метод второго порядка. Согласно этому методу, на каждом итерационном слое вместо функционала (17) применялся функционал
JN (Л+1) = J (л) + и(л+1 — Л')
(25)
где N - достаточно большая положительная постоянная.
На рис. 1 показано поведение определяемого параметра в зависимости от величины N . При этом уравнения (23) и (24) решались с помощью неявной разностной схемы. Сетка разбивала пространственный отрезок х е [0,1] на 10 интервалов, временной отрезок t е [0, 0.8) на 100 интервалов. Точное решение вычислялось в 5 х 50 точках пространства-времени. Кривые 1 -4 соответствуют следующим значениям множителя
х
х
4
1
Рис. 1. Сходимость определяемого параметра: 1 - N = 100 ; 2 - N = 40 ; 3 - N = 20 ; 4 - N = 10
N = 100; 40; 20; 10 . Как видно, введение модифицированного функционала (25) очень напоминает введение демпфера в колебательный процесс. При достаточно больших значениях «вязкости» (N = 100) процесс перестает быть колебательным и параметр X подходит к точке равновесия с одной стороны. При малых значениях N (Ы = 10) параметр X совершает незатухающие колебания вокруг точки равновесия.
Влияние количества «измерений» на качество процесса восстановления параметра X показано на рис. 2. В этом случае уравнения (23) и (24) решались на сетке 10 х 100 узлов при N = 40. Кривые 1-4 соответствуют случаям, когда точное решение вычисляется в 5 х 50, 3х 50, 5 х 20, 10 х 50 точках пространства-времени. Как видно, количество пространственных измерений довольно существенно отражается на качестве восстановления. Увеличение же количества замеров по времени ведет к повышению чувствительности, процесс восстановления приобретает колебательный характер.
В связи с тем, что при решении уравнений (23) и (24) используется разностная аппроксимация, полученное в результате значение искомого параметра несколько отличается от «точного» значения X = 0,2 . Таким образом, фактически имитируется случайная погрешность наблюдений. Вообще говоря, такая ситуация имеет место всегда, так как при использовании конечно-разностных методов идентифицируются не параметры исходного уравнения, а параметры его разностного аналога.
Рис. 3 демонстрирует влияние количества временных слоев разностной схемы на качество восстановления. Точное решение вычислялось в 5 х 20 точках пространства-времени при N = 40. Кривая 1 соответствует сетке из
Рис. 2. Влияние количества замеров на качество сходимости итерационного процесса: 1 - 50^50 точек; 2 - 3^50 точек; 3 -5x20 точек; 4 - 10x50 точек пространства-времени
(26)
Рис. 3. Влияние количества временных слоев разностной схемы на сходимость итерационного процесса: 1 - 100 временных слоев; 2 - 40 временных слоев
10 х100 узлов, кривая 2 - сетке из 10 х 40 узлов. Как ни странно, сетка с меньшим количеством временных слоев (кривая 2) приводит к более точному восстановлению параметра X. Это, по-видимому, связано с накоплением погрешностей при увеличении количества вычислений. Заметим, что применение сеток из 10 х120 и 10 х 20 узлов дает кривые восстановления, практически совпадающие с линиями 1 и 2, соответственно.
Теоретически необходимость введения модифицированного функционала (25) может быть обоснована с помощью теорем теории оптимального управления системами, описываемыми уравнениями с частными производными параболического типа, сформулированных Ж.-Л. Лионсом [4].
Уравнение (22) на верхнем итерационном слое может быть переписано в виде
= х д 2и*+1
дt
дх2
д 2и
+ f, (26)
дх
а функционал (25) в виде
JN (дА)=^^\[ы(х],^)-г, + N (ЗА)2. (27)
]=1 о
Существует единственное значение (ЗА), удовлетворяющее условию минимума функционала (27) на (5 +1) -м итерационном слое только при условии N > 0 [4]. Если N = 0, то условие единственности может быть не выполнено. Для того чтобы условие единственности было выполнено и при N = 0, требуются измерения производных по времени от функции и . В этом случае функционал задачи должен иметь вид
тЫи(х],%) (Л 2
jn=j+Zi
j-i 0
dt
d%, (28)
где z1 - дополнительные данные наблюде-
нии.
Однако на практике непосредственное измерение производных затруднительно, а дифференцирование данных наблюдений функции и вносит большие погрешности.
Следует отметить, что условия минимума модифицированного функционала (25) совпадают с условиями для исходного функционала (17). Интересно, что введение в исходный функцио-
нал (17) добавки вида (25) аналогично введению своеобразной функции штрафа на каждой итерации. Мы как бы «штрафуем» функционал за большое отклонение нового значения As+1 от значения As на предыдущей итерации.
Как показывает опыт решения конкретных задач, привлечение дополнительной информации о решении (знание производных или более разнообразный набор измеряемых функций в случае систем уравнений) улучшает качество процесса восстановления неизвестных параметров.
Литература
1. Алифанов О.М. Определение тепловых нагрузок из решения нелинейной обратной задачи // Теплофизика высоких температур. - 1977. - Т. 15, № 3. -С. 598-605.
2. Chayent G., Dupuy M., Lemonnier P. History matching by use of optimal theory // Soc. Pet. Eng. J. -1975. - V. 15, № 1. - P. 74-86.
3. Seinfeld J.H. Identification of parameters in partial differential equations // Chem. Eng. Sci. - 1969. - V. 24, №. 1.
4. Лионе Ж.-Л. Оптимальное управление системами, описываемыми уравнениями с частными производными. - М.: Мир, 1972. - 414 с.
Поступила в редакцию 23.10.2014
