Кинематика фильтрационных течений и ее приложение для решения обратных задач гидропрослушивания Текст научной статьи по специальности «Математика»

УДК 550.833

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

А.И. КОБРУНОВ

Ухтинский государственный технический университет, г. Ухта [email protected]

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

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

A.I.KOBRUNOV. KINEMATICS OF FILTRATION FLOWS AND ITS APPLICATION FOR THE SOLUTION OF INVERSE PROBLEMS OF WELL INTERFERENCE TESTING

We considered the kinematics of filtration equation as the regularity of characteristic point movement of pressure restoration curve. The equation is written in the form of integral on trajectories for interval times of backoff point distribution of the pressure restoration curve between injection wells and reception wells. On the basis of the constructed equation we carried out the statement of an inverse task for interference testing to find the spatial distribution of piezoconductivity coefficient within the permeable seam in the form of tomographic task.

Key words: well interference testing, filtration flows, piezoconductivity coefficient, inverse tasks, tomography, optimization theory, program and algorithmic software, risk decrease

Введение

Пусть система из N скважин, вскрывших гидродинамически связанные между собой участки проницаемого пласта, допускает проведение следующего эксперимента. Выберем некоторую скважину из N , скажем с номером і, подвергнем ее длительной депрессии (положительной или отрицательной) в пределах продуктивного пласта. В остальных скважинах или части из них числом М < N будем регистрировать время г , і = 1 + М , за кото-

і ^ ^ і

рое депрессия дошла до них и достигла максимальной скорости - точки перегиба нарастания (или спада) давления. В итоге будет получено Мі

значений времен для веера распространения точек перегиба на кривых изменения давления от скважины і к остальным скважинам. Поменяем точку і на і +1 и повторим эксперимент. Получим новый веер значений г1+1], і = 1 ^Мм . Продолжая процесс, получаем систему значений интервалов движения экстремальной точки (точки перегиба) между

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

Дифференциальные уравнения

Параметрами, входящими в формулировку уравнений подземной гидродинамики, служат: ско-

(1)

рость фильтрации v(t, x) = {vx, vy, vz}; плотность

заполняющего флюида р(t, x), x = {x, y, z}; t -

время, давление p (t, x) и пористость m (t, x). Исходным уравнением для перечисленных компонент служит хорошо известное уравнение неразрывности [2]:

, , ч , чЧ d[p( t, x) m (t, x)]

div (р (t, x) v (t, x)) = - [^ ^ y n

Оно справедливо в линейном приближении к процессам движения флюидов в отсутствии внешних источников. Если источники присутствуют и их

величина равна q(t, x) , то они добавляются как

слагаемое в правую часть уравнения (1).

Основным постулатом в теории фильтрации служит принятие закона о скорости фильтрации в зависимости от градиента давления

grad (p (t, x)) и некоторой массовой силы f : grad (p (t, x )) + f = ® (| vI) v / |v| . Функция

ф (I v|)зависит от модуля скорости фильтрации и

определяет силу сопротивления при движении в пористой среде. Если эта функция линейна, напри-

мер

ф (I vl)

v , где: /и - динамическая вяз-

кость, k - коэффициент проницаемости пористой среды, то закон называется линейным законом Дарси, а текущая вязкая жидкость - Ньютоновской жидкостью. Нетрудно видеть, что закон связи скорости фильтрации и градиента давления p (t, x) в этом, последнем случае ньютоновских жидкостей,

примет вид: v = - [ grad ( p (t, x )) + f

Роль

V = -

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

Фундаментальное свойство Ньютоновской жидкости состоит в независимости ее вязкости и

как следствие силы сопротивления Ф V]) от градиента скорости движения. В противном случае жидкость неньютоновская. Подставляя полученное выражение для скорости в закон сохранения, получим для ньютоновской жидкости уравнение движения:

div f р(t,x) и [ grad (p (t,x ))]^=-5[p( t ’(t ’x ^

(2)

Величина K =— называется динамическим

и

коэффициентом фильтрации,

L3T

имеет размерность , в то время как коэффициент проницаемости

М к: [ С ].

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

Плотность жидкости и коэффициент пористости зависят от давления через коэффициенты сжимаемости жидкости (Ж и скелета (С :

Р

Р

dm = PCdi

(3)

или:

grad р (t, x) = РЖ р (t, x) gradp (t, x);

Dm _ dp

~dt ~ ^ ~dt'

(4)

Объединяет коэффициенты сжимаемости жидкости (Ж и скелета (С коэффициент упруго-

емкости пласта ( = т(Ж + (С , который определяет изменение упругого запаса жидкости dVЖ в объеме V: dVЖ = (*У0dp .

Коэффициент пьезопроводности К выражается через коэффициенты динамической фильтрации К и упругоемкости по правилу

K

k

к = -

T

(5)

Пользуясь коэффициентом пьезопроводности и проведя ряд преобразований [3], получим:

div [к (x) (gradp(t, x))] = p( , ) , ему эквивалентное

div [к (x) (gradp (t, x))] = p ( , ) .

(б)

(7)

L \_p(t, x )] = f (t, x);

p( t, x L=g(t, x);

p( t, x )l t=o=Po(x).

(її)

Эквивалентность этих уравнений определена линейной связью между плотностью и давлением. Они служат основой для последующего анализа.

Приведенные уравнения характеризуют «свободное распространение» без каких-либо внешних побуждающих факторов. Чтобы это движение было не тождественным (состоянием покоя), должны присутствовать некоторые начальные условия р0 (x) , характеризующие состояние системы в некоторый выбранный момент или моменты времени в области D , краевые условия g (t, x) , характеризующие режим распространения давления на границе dD той области, где это распространение происходит и, наконец, источники или стоки f (t, x) , действующие мгновенно или на протяжении некоторых интервалов времени внутри D .

1.Функции Грина. Для построения решений уравнений (6) и (7), что с формально математической точки зрения одно и то же, необходимо дополнить их краевыми и начальными условиями, а также сделать предположения о поведении коэффициента пьезопроводности. Первое такое предположение состоит в том, что к(x) = const. и уравнения (6)-(7) трансформируются в:

",x ^ (8)

к дt

В цилиндрической системе координат для

р( r, ©, z ):

д2 1 д 1 д2 д2

2 +--------+—2------------------------2 +-2 (9)

дr r Sr r д© дz

r = -у/ x2 + y2; © = arctg У; z; x = r cos (©); y = r sin (©)

В сферической системе координат для

P(r, ©,х):

, x )|t=o = р0 (x)

Существует стандартизующая обобщенная функция w(t, x) , линейно зависящая от левых частей в (11) такая, что (11) эквивалентная задаче [4]:

L\p(t, x )] = w (t, x );

p(t, x )Ld=(; (12)

P(t,x )l t=o = 0

Для стационарных задач, уравнения которых инвариантны относительно сдвига во времени, а уравнение (б, 7) относится именно к таким задачам, решение (12) задается соотношением:

t

P(t, x ) = ff G (t-т,x -\)drd^ . (13)

0 D

Функция G называется функцией Грина для уравнения L [G (t, x)] и находится как решение задачи:

L\G (t, x )] = 8( t )8( x );

G (t, x L=(; (14)

G (t, x)lt =o = (.

Здесь 8 ( ) - дельта функция Дирака.

Для уравнения

др(^ x)

- кДр(t, x ) = f (t, x);

дt

P(t, x )l t=o = Po (x ) ,

(15)

дадада

|p(t,x)|dx < C

-да -да -да

д2 , 2 д i 1 д2 + ctg х д

дr2 r дr r2 дх2 r2 дх

(10)

И С не зависит от ^, в бесконечной области для t > 0 и начальном условии <р(7,= (р0 (х) стандартизующей функцией служит [4]:

ы (t, х) = / (t, х) + р0 (x)J(t) . (16)

Функция Грина:

П-----2----2 ъ У \x" + У

r = у]x + У + z ;© = arctg —;х = arctg—----------

G (t, x ) = -

1

exp

2\J кжt

x = r cos (0) sin (x); y = r sin (0) sin (x); z = r coc (x).

Для дифференциального уравнения, дополненного действующими источниками, начальными и краевыми условиями

(17)

Здесь как обычно, г = х = х+ у + г и

х - \ = і(х - ) + і(у - ) + k(г - ).

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

з

ния. В том случае, если источник работал некоторое время - 5 с постоянным дебитом, равным единице, распределение давления в пространстве вычисляется в соответствии с (13):

s exp

р( s, t, x ) =

4к(і -т)

____1____f

(V4кж) o (t-т)

4т. (18)

В том случае, если движение осуществляется в горизонтальном пласте, на кровле (г = 0) и подошве (2 = I) которого поддерживаются условия р(1,х)|г=о = & (х,у,t);р(t,х)|г=1 = g2 (х,у,t), то задача (15) трансформируется в

др( t, x )

ді

- кДр( t, x ) = f (t, x );

р О1, х )|= Р0 (х ), (19)

р(t, х )1 г=0=&1 (х, y, 1:) ;р(1 , х )1 г==&2 (- y, 1:).

Стандартизующей функцией служит:

ы ^ х ) = I ^, х) + Р0 (х )8^)-к8\г) & (х, у, t )-~К8' (I- г) &2 (х, y, t),

(20)

где 8'(г) производная от дельта функции Дирака,

а

определяемая условием: | ^(г8'(г)dz = (0), а

-а

ее решение определено функцией Грина:

G ( ^ У, z,^x ,4y ,4z, i):

2

l (2>/ кжі )3

exp

(x -£x )2 + (У -4, )2 4кі

. Knz . жп4 Г кп1ж1і

x> sin---sin----- exp I----—

S і і 4 і ,

(21)

2. Пьезокинематика. Для прямолинейно параллельного движения, когда единственной существенной координатой является координата x , уравнение (8) трансформируется в

д2 ч 1 др( t, x)

— р(t, x ) = -

x ' ' к ді

а функция Грина для него равна:

1

р(t, x)—pexp Vi

x2

4кі

(22)

(23)

Для плоско радиального осесимметричного относительно вертикальной координаты движения флюида, которое соответствует выбору цилиндрической системы координат и в которой движение симметрично относительно оси OZ и не зависит от угловой координаты 0, уравнение (8) трансформируется в:

д

1д

^ p(t,r ) + ZT: р( t,r ) = 7

1 др(і, r)

дr

r дr

ді

(24)

Функция Грина для этого уравнения есть:

p(t, r) =1 exp

4кі

Наконец, для сферически радиального дви-

жения:

д

2д

^ p(t,r p(t,r ) = 7

1 др(і, r)

дr r дr

Функция Грина:

к ді

p( t, r) = —= exp t\t

4кі

(25)

(2б)

Функция р(t, г) во всех трех случаях (23, 25,

27), где в первом случае г = х, имеет одну и ту же точку перегиба как функция переменной г . Действительно, все эти уравнения, следуя Щелкачеву [3], можно записать в единой форме:

P( t, r ) =

1

1+а

t ~Т

exp

r2 4кі

(27)

Случай (23) соответствует а = 0, (25)-

а = 1 ,и, наконец (27)-а = 2 . Дифференцируя последнее выражение дважды, по г получим:

д2р(і, r )_ 1

f 2 Л 3+а

дr2

2к

2кі

-1

t

exp

f \ 4кі

= 0.

/

Откуда уравнение для экстремальной точки

Г ті

2кі

-1

= 0 и,следовательно:

Ге =л/2К1 . (28)

Это уравнение для точки перегиба. Его следует интерпретировать следующим образом. Для каждого момента времени tкривая давления р

имеет точку перегиба на расстоянии г = V2к1 от

начала координат. Иными словами, местоположение точки перегиба есть функция времени. Следовательно, в координатах (ге, I) точка перегиба

движется в пространстве и скорость этого движения

V St 2к к e = ді _ 2у[2ш ~ re

(29)

Точка перегиба для пространственного распределения функции р(1, г) при зданном времени

t фиксирует пару координат |ге,!} . Для выбранного

момента времени, соответствующее ге - это точка

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

время достижения максимума давления, после которого оно начинает спадать. Для того, чтобы убедиться в этом, найдем максимум (27) по t при

2

r

2

r

2

2

r

2

x

2

фиксированной точке г. Для этого продифференцируем (27) по времени и приравняем результат нулю:

др(і, r ) ді _

Но тогда

r2

1 + а

Л

4кі 2

3+а

2

2

exp

4кі

= 0.

і = , (30)

2к (1 + а)

что совпадает при а = 0 с полученным выше выражением (28) для координаты ге точки перегиба. В

данном случае, это уже пространственная координата, в которой через найденное время і будет наблюдаться максимальное давление. Рассматривая скорость Vm (г) перемещения точки максимума давления в пространстве, получим:

дг ду] 2кі (1 + а) 2к (1 + а) к (1 + а)

V (r) = —=

m

ді

ді

2^2кі (1 + а)

(31)

При а = 0, что соответствует одномерному

плоскопараллельному движению, Vm (г) = Ve (г).

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

Приведенные соотношения относились к случаю мгновенного источника, включенного в начале координат. Если перейти к источнику с постоянным дебитом, расположенным в начале координат и действующем на протяжении времени і, это соответствует вычислению интеграла (18) для 5 = і . Рассматривая фундаментальные решения для плоскопараллельного, плоскорадиального и сферически симметричного потоков, получим:

P(t, r ) = f

1

1+а (i-р-

exp

r

4к (t-т)

(З2)

Вычисляя вторую производную от р(і, г) по

і, и приравнивая результат нулю, получим уравне-

„2

ние для точки перегиба [3, с. 140]: і =

2к (1 + а)

Это в точности выражение (30), но уже не для экстремума от точечно и мгновенно действующего источника, а точка перегиба кривой восстановления давления от источника, представляющего собой постоянный в течение времени t сток или приток. Повторяя приведенные выше рассуждения и нахо-

дя скорость движения точки перегиба кривой восстановления давления, получаем ту же формулу (31)

кп

V ( r ) = —

(ЗЗ)

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

Найденное выражение (33) для скорости распространения характерных точек восстановления давления, позволяет перейти к понятию пространственно распределенного поля времен т(ro, r)

распространения характерных точек восстановления давления регистрируемых в точках r и вызванного точечным источником расположенном в r. Выберем точку r' и построим проходящую через нее поверхность w(t,г) как изохрону для

т(r') = const. Если через точку r' изохрона проходит, то она только одна. Любое другое движение приходит в эту точку за большее время, либо не приходит вовсе. За интервал времени At расстояние, которое пройдено характерной точкой давления от jr'} , образует поверхность сферы радиуса

Ar = V ( r) At;

r + Ar

At =-------.

кп

Огибающая этих сфер образует новое положение изохроны с моментом времени T + At . Продолжая этот процесс и так далее, получаем для времени движения от точки r0 до произвольной точки r правило:

l (£) d$

(З4)

I(%) представляет собой участок кривой С(г0,г) от точки г0 до текущей точки % .

Линия С(г0,г) минимизирует (34) по всем

возможным кривым, соединяющим г0 и г. По этой

причине она относится к стоящему в знаменателе под знаком интеграла выражению коэффициента

пьезопроводности и обозначается также Ьк (г0, г) .

Следует обратить внимание на то обстоятельство, что в отличие от уравнений полей времен для гиперболических уравнений в (34) - для параболических задач, время движения от промежуточ-

t

2

r

ной точки до конечной зависит от траектории, по которой пришло возмущение в промежуточную точку. Это обстоятельство накладывает ограничения на класс алгоритмов для поисков траектории минимизирующей (34).

3. Постановка обратной задачи гидропрослушивания. Впервые задача нахождения коэффициента пьезопроводности пласта по результатам анализа кривой восстановления давления между двумя скважинами по сведениям Щелкачева [5, с. 256] была сформулирована Б.С. Харченко в 1957 г., и далее развивалась Г.И. Баренблатом, В.М. Енотовым и В.М. Рыжиком в 1984 г. В этих работах развивалась методика обработки измерений в регистрирующей скважине под влиянием нагрузок в закачивающей скважине. Определялось одно числовое значение, характеризующее величину коэффициента пьезопроводности для всего участка между скважинами. Попытки нахождения пространственного распределения этого коэффициента не делалось. Уравнение (34), примененное для расчета серии интервальных времен движения характерных особенностей кривых пластовых давлений по сети из серии скважин источников - приемников, позволяет приступить к постановке обратной задачи гидродинамического прослушивания с целью нахождения пространственного распределения коэффициента пьезопроводности.

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

площади входящей в пределы участка залежи с гидросвязанными между собой элементами будем считать, что заданы результаты реальных измерений в парах скважин числом т = ^, ]) по системе

точек , г } :тг (г, г]) . Они соответствуют некото-

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

к0(%). Пусть расчетное время распространения

особенностей кривой восстановления давления по (34) для заданного коэффициента пьезопроводности

.(%)

есть Т

(г, г. ) = тп (г, г. ).

Ч)\ ' ^ ) 0\ ^j/

Обозначим Дк0 (% ) = кг (%) — к0 (% ) и Дт0 (г, г ) = т0 (г, г ) — тг (г, г ) . Примем, что величина Дк (%) имеет малые вариации и приближенно можно считать Ьк (г0,г) = Ьк (г0,г) . В этом

случае:

Дт (гг)- Г 1к0%%% — Г !лШ

0 (:''Ь Ч(I,I3к. (%) Ц,)*' %

'., (%)

, к (%)— г

Г 3к (%) Г

Ч (г г

3к0 (%) ^^] 3[к0 (%) + Дк (%)]

'„ (%)[Дкс (%)]4%

3[к0 (%) +Дк0 (%)]хк0 (%)

'0 (%)[Дк0 (%)]Л%

3к02 (%) '

Таким образом, для невязки времен прихода возбуждения получаем линеаризованное для (34) уравнение:

Дт0 (г.г, ) = т0 (г. г )-Т (г.г,)

Г

'„, (%Ж (%)]

3к2(%)

(35)

Уравнение (35) - это линейное интегральное уравнение относительно добавки Дк0(%) к нулевому приближению для пространственного распределения коэффициента пьезопроводности. Его решение осуществляется по стандартной схеме сведения линейных интегральных уравнений к системам линейных алгебраических уравнений либо на основе метода конечных элементов, либо его частного случая - метода сеток [6]. В рамках выбранной параметризации для Дк0 (%) , вектор искомых

параметров Лк0 связан с вектором невязок

Ат0времен прихода возбуждения матрицей А,

рассчитываемой на основании интегрирования по траекториям (35) уравнением:

ААк0 = Ат0 (36)

Таким образом, задача (35) сводится к (36), и эта суть - основа для построения приближенных устойчивых решений обратной задачи [7, 8] гидродинамического прослушивания. В процессе ее решения должны быть использованы итерационные схемы уточнения решений и методы.

Г

Г

к

Литература

1. Терещенко СА. Методы вычислительной томографии. М.: ФизМатЛит, 2004. 319 с.

2. Басниев К.С., Дмитриев Н.М., Каневская РД, Максимов В.М. Подземная гидромеханика. М.; Ижевск: Институт компьютерных исследований, 2006. 487 с.

3. Щелкачев В.Н. Основы и приложения теории неустановившейся фильтрации. Часть 1. М.: Нефть и газ, 1995. 586 с.

4. Бутковский А.Г. Характеристика систем с распределенными параметрами. М.: Наука, 1979. 219 с.

5. Щелкачев В.Н. Основы и приложения теории неустановившейся фильтрации. Часть 2. М.: Нефть и газ, 1995. 586 с.

6. Марчук Г.И. Методы вычислительной математики. М.: Наука, 1980. 534 с.

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

8. Кобрунов А.И. Математические основы теории интерпретации геофизических данных. М.: ЦентроЛитНефтеГаз, 2008. 286 с.