CC BY
32
УДК 517.958:532.546+556.332.52+556.334+517.551+517.584
АНАЛИТИЧЕСКОЕ ОПИСАНИЕ УСТАНОВИВШЕГОСЯ ВОЛНОВОГО ГИДРОГЕОДИНАМИЧЕСКОГО ПРОЦЕССА В НАПОРНОМ ПЛАСТЕ С ПЕРЕТЕКАНИЕМ
Т.Ю. Заведий, Л.Я. Ерофеев*
ОАО «Сибирский химический комбинат», г. Северск *Томский политехнический университет E-mail: taras_zavedy@mail333.com
Выполнена практическая аппроксимация модифицированной функции Бесселя второго рода нулевого порядка (функция Макдональда) и производных от нее расширенных функций комплексного аргумента применительно для случая решения задачи волновой гидрогеодинамики об установившихся гармонических колебаниях уровня подземных вод в бесконечном по простиранию напорном пласте с перетеканием. Показана принципиальная возможность упрощения сложного решения с перетеканием и сведения его к решению более простой аналогичной задачи для изолированного пласта.
Ключевые слова:
Волновая гидрогеодинамика, установившийся колебательный режим, напорный пласт с перетеканием, асимптотические аппроксимации функций Бесселя, функция Макдональда от комплексного аргумента.
Key words:
Wave hydrogeodynamics, steady oscillatory regime in aquifers, confined aquifer with leakage, asymptotical approximations of Bessel functions, McDonald function with complex input values.
Бесконечный по простиранию напорный пласт с перетеканием с точки зрения аналитической теории решения прикладных гидрогеологических задач является, пожалуй, одним из сложных случаев [1, 2], причем дальнейшее усложнение гидрогеоди-намической модели такого пласта зачастую требует уже применения численного моделирования на ПЭВМ в специализированных приложениях. С этих позиций, в качестве цели работы, интересно рассмотреть получение аналитического решения и практическое упрощение этого решения для задачи о затухании амплитуды установившихся гармонических колебаний, возбуждаемых на стенке скважины [3] в бесконечном напорном пласте постоянной толщины при наличии перетекания.
Дифференциальное уравнение для приращения напора подземных вод £(г,/) при таких условиях примет следующий вид [1, 2]:
1 as (r, t)
a dt
(
d2S(r, t) 1 dS(r, t)
dr2
-+ — r
dr
- b2S(r, t), (1)
Подстановка (2) в (1) после сокращений приводит для функции затухания комплексной амплитуды А(г) к уравнению эллиптического типа [4]:
d2 A(r) + 1 dA( r)
dr
m
+ b2 I A(r) = 0,
(3)
с граничными условиями на стенке скважины и на бесконечности: А(г)г=гл=А1, А(г)м„=0.
Дифференциальное уравнение (3) путем простых преобразований сводится к удобному безразмерному виду:
1 ЙА
Р йР
d2 A
~Т + _“—(i + Р)A = 0 dp
(4)
где p = rjrp = r.— - безразмерный радиус от
ис-
где г - радиус от источника волнового возмущения (скважина), м; / - время, с; г0 - радиус стенки скважины, м; А1 - начальная амплитуда колебаний в
, % 1
скважине, м; ь = л^ - параметр перетекания,
1/м, обратно пропорциональный фактору перетекания В; % - коэффициент перетока разделяющих водоупоров, 1/с; Т - водопроводимость пласта, м2/с; а - пьезопроводность пласта, м2/с; со - угловая частота периодического процесса, с-1. Граничные условия на стенке скважины и на бесконечности: 8(г,Р)г==А1$,т.Ш, £(г,/)м„=0.
Решение (1) для установившегося гармонического процесса без начальных условий ищется в комплексном виде [3, 4]:
5 (г, *) = А(г)ёш‘. (2)
точника волнового возмущения; гр - характеристический радиус влияния периодического процесса; 3 - волновой параметр перетекания, имеющий простую интерпретацию - это квадрат отношения Р=(гр/Б)2характеристического радиуса влияния периодического процесса и фактора перетекания.
Решение (4) известно [4, 5] и записывается в виде:
А(Р) = С^(4Т^-р) + С210ЦI + Р-р), (5)
где 10(х) - модифицированная функция Бесселя первого рода нулевого порядка (функция Инфель-да); К,(х) - модифицированная функция Бесселя второго рода нулевого порядка (функция Макдональда); С1, С2 - постоянные интегрирования, из граничных условий (1) и асимптотических свойств 10(х) очевидно, что С2=0.
В математической физике известно, что функ-ц-и-я Макдональда от комплексного аргумента вида ЛХ может быть представлена двумя отдельными функциями от вещественного аргумента х, так называемыми функциями Томсона (Кельвина) -кег(х) и кеЦх) [4]:
К0(л/Ї • х) = кег(х) + і • кеі( х). (6)
Тогда, по аналогии с (6) для вычисления и анализа полученного решения (5), необходимо ввести расширенные функции Томсона:
К0(7і + Р • х) = Кег(х,Р) + і• Кеі(х,Р) (7)
причем очевидно, что при р=0 расширенное выражение (7) переходит в (6).
Окончательно решение (1), с учетом (5) и (7), для варианта экспоненциального представления комплексного числа запишется в виде:
I Кег (р, Р) + Кеі (р, Р)
КеГ2(Ро, Р) + Кеі2(Ро, Р)
Кеі(р,Р)
КегрР
Вег(х, Р) |^С+1п^^ + Р2х І I-
-Веі( х, -Р
1 ^ 2 J
Х тХ (-1)кСк Рп
(9)
Кеі( х, Р) = Зш[К0(^/ і + Рх)] = Веі( х, Р) ^С + 1п [ ±^х Ц +
+Вег( х, Р)-Таге1§ -Р
Х
1
-(2к +1)
(10)
где скобки [и/2] обозначают операцию округления и/2 до ближайшего меньшего целого;
Вег( х, Р) = ^е[10(л/^"+Рх)] =
=Х
гХ (-1)кС2к Рп
'= 22 п («!)2^
Веі( х, Р) = Зш[І0^ і + Рх)] =
= Х:
Х (-1)ксПк+іРп
-(2к+1)
хе [. {^Ю)}. (8)
где р0=г0/гр - безразмерный радиус стенки скважины.
Собственно, в дальнейшем, интерес представляет только первая подкоренная часть выражения (8), поскольку именно она характеризует максимально возможную амплитуду колебаний, а комплексная экспонента - вторая часть выражения (8), в таком случае, для решения задачи (1) представляет меньший интерес, поскольку описывает пространственное запаздывание фазы колебаний и сами периодические колебания во времени.
Если с вышеприведенными выкладками по получению аналитического выражения решения (1) не возникло значимых затруднений, то практический расчет расширенных цилиндрических функций Томсона представляет довольно серьезную задачу вычислительной математики. Чтобы проиллюстрировать всю сложность, достаточно привести выражения рядов Тейлора-Маклорена для введенных в этой работе расширенных функций Томсона, которые легко выводятся с учетом аналогичных выражений для кег(х) и кеі(х) из [5, 6]:
Кег( х, Р) = 0( і + Рх)] =
1
и=0- п(и!) *=0
Из общей специфики вычисления значений функции Макдональда следует отметить, что ряды (9, 10) для задач прямых математических расчетов практически бесполезны, поскольку требуют учета очень большого количества членов ряда и фактически обеспечивают приемлемую сходимость для значений аргумента х не более 4...4,5. В подобной ситуации необходимо использовать для получения новых вычислительных аппроксимирующих выражений асимптотические свойства цилиндрических функций, что широко применяется в различных расчетных схемах [7].
Но с другой стороны, прежде всего для задач геологических дисциплин, часто возникает необходимость получения еще и максимально простых выражений аппроксимации специальных функций математической физики, что обычно требуется для интерпретации и оценочных полевых расчетов, производимых с невысокой точностью, без применения сложных вычислительных средств и алгоритмов. Одним из таких примеров компромисса между точностью и простотой задания аппроксимации является полученное авторами работы базовое кусочно-аппроксимирующее выражение цилиндрической функции Макдональда К0(х) для двух интервалов задания значений аргумента (0,1), [1,<х>] с простым дробно-рациональным заданием коэффициентов при степенях аргумента и удобным максимально точным сопряжением в окрестности аргумента х=1.
Ко(х) *
-|С + 1п2
-х Ґ
1+
4~х
1 —
117 107 ^
1
2
х< 1;
V
х > 1.
(11)
(8х + 7/2),
Модуль абсолютной погрешности сопряжения в точке х=1 равен 1,19761-10-5. Значения функции, рассчитанные по (11), сравнивались сточными значениями для К0(х) из [8]. Результаты сравнения приведены на рис. 1.
Из рис. 1 видно, что максимальная относительная погрешность аппроксимации не превышает 1,45 %, для интерпретационных решений гидрогеологических задач это вполне приемлемо. Максимум погрешности аппроксимации лежит в интервале значений аргумента х [0,5.0,9]. Это об-
Рис. 1. Зависимость ошибки аппроксимации К0(х) отаргумента х
условлено тем, что уменьшение погрешности аппроксимации (по асимптотическим свойствам) наблюдается у аппроксимированной функции Ко(г) при х^0 либо при х^да. За исключением интервала значений аргумента [0,17.1,0], модуль относительной погрешности аппроксимации не превышает 0,1 %, а в указанном интервале относительная погрешность аппроксимации не превышает 1,45 %.
Далее, из аппроксимирующего выражения (11), можно получить комплексное решение для вещественной и мнимой частей К^) при произвольном комплексном значении аргумента 1=У+т для первого выражения из (11) при (0<|^|<1):
Ко(2) *-\С + 1п-
2 Л
1 + — £_ 107 4
1 +117 (V2 - м2 + 2г'гм) | + 428 )
2 2 -.ч
+------(V - м + 2^и>) =
428
С + 1п2л/V2 + м2 |>
-г
л 117 , 2 2, 1 117 Л м
1 +---------(V - м ) I----2vw аг^ —
428 I 428 V,
117(2 2)
-----(V - м ) -
428
2ум ^С + 1п-1-у/V2 + м2 | +
(л 117, 2 ^ 1 м 117 „
+ 1 1 +-(V - м ) I аг^-----------------2vw
I 428 I V 428
(12)
Аналогичным образом получено комплексное решение для вещественной и мнимой частей К0(г) из аппроксимационного выражения (11) при произвольном комплексном значении аргумента 1=у+т для второй части выражения при |^|>1.
К0( 2) ^\1 --
1
2 42 \ (82 + 7/2)
п е2 ( 16 г + 5
2 42 \ 16 2 + 7
-)
2 4(у + гм)
(256(v2 + м2) +192v + 35) + г'32 м
Ке[К0(2)] *
-(256^2 + м2) + 192v + 35) х
( 1 м 1
хсоб1 м + — arctg — 1 +
I 2 VI
I 1 м
+32м• Б1п| м + ^гС^ —
• е
((16v + 7)2 + 256м2)4у2+м2,
Зш[К0(2)] *
^'-(256^2 + м2) +192v + 35) х
• ( 1 + м 1
хБ1пI м + — arctg — 1 +
I 2 V I
((16 V + 7)2 + 256
(13)
(14)
Из (13) получены выражения (14), (15) для мнимой и вещественной частей. В дальнейшем удобно будет видоизменить аргумент расширенных функций Томсона в несколько иной форме записи для входных переменных, чем в (9), (10): Ке1*(у,^)=Зш[Кв(г)], Кег*(у,^)=ЭТе[Кв(7)], тогда амплитудная часть решения (8) запишется в следующем безразмерном виде:
А(р) = \К р(2 р)| =
А К, (20, Р)|
( Ker*2(v(p, Р), w(р,Р)) +Л +Kei*2(v(р, Р), w( р, Р)) Ker*2(V(p0, Р), w(p0,Р)) + +Kei*2(v(po, Р), w(ро, Р))
(16)
где
z = v + iw
= yfi+Pp =
p^1 + Р2 cosfiarcctg^) J +
+1р^1 + Р2 sin^-2агсс1§(Р)J;
zo = vo + lwo =41+ M =
( 1
и в случае р=0, только с большим наклоном на логарифмическом графике. Т е. вычисление корня из суммы квадратов расширенных функций Томсона (16), по крайне мере при 0<Р<1, всегда можно аппроксимировать более удобными для расчета функциями Томсона - кег(к(Р)-р) и ке1(к(Р)-р) с коэффициентом масштабирования к(Р) входной переменной р, что в несколько иной форме записи будет выглядеть как:
\к0(ру1Г+Р)\ * \к0(к(Ррг)|. (17)
Графики решения (16) с использованием аппроксимации (12, 14, 15) расширенных функций Томсона при различных задаваемых значениях волнового параметра перетекания Р изображены на рис. 2. Т е. результаты расчета (16) довольно легко были аппроксимированы более простым выражением (17) с максимально точным совмещением кривых на графике с логарифмическим отображением оси ординат. Другими словами, аппроксимирующее упрощение решения (1) для определенного авторами работы выражения корректирующего коэффициента масштабирования к(Р) аргумента с учетом (16, 17) будет выглядеть как
А(р) = \Ko(P^^i + p) |K„(p-(1 + 0,5158Р))|
А
ро\А + Р2 cos(^arcctg(Р) J + +ро^РТ-^ ^rcrtg^) J; v(p, Р) = рф + Р2 cos ^^2arcctg(Р)J;
w(p,Р) = рф + Р2 sin^^arcctg^)J.
Ниже приведен пример расчета затухания гармонических колебаний, возбуждаемых в одной скважине в зависимости от увеличения радиального удаления. Для типичных значений радиуса многих обсаженных скважин, равных в среднем 0,1 м и для значений характеристического радиуса периодического процесса, например равного rp=7,98 км (при периоде колебаний 400 суток и пьезопроводности я=1-106м2/сут), безразмерный радиус стенки скважины p=r0/rp будет равен 1,25-10-5. Поэтому в знаменателе выражения (16) уверенно можно использовать аппроксимирующие выражения для модуля аргумента z меньше 1 (12).
Результаты расчета по (16) для интервала задаваемых значений волнового параметра перетекания 0<Р<1 дополнительно выявили одну очень важную особенность у полученных графиков. С увеличением значения волнового параметра перетекания Р радиальное затухание колебаний происходит быстрее, но самое главное, по такой же подобной функциональной зависимости, как
K^^) |K„(p„ -(1 + 0,5158Р))|
при 0 < Р < 1.
На рис. 2 в иллюстративных целях специально показан рост погрешности сопряжения в точке р=1 при увеличении значения волнового параметра перетекания Р. В дальнейшем, при разработке вычислительных алгоритмов аппроксимации функции (7), необходимо всегда предусматривать коррекцию коэффициентов аппроксимирующего выражения для лучшего сопряжения двух решений в точке р=1. Это особенно актуально при больших значениях Р>1. Однако эти вопросы необходимо будет рассматривать в отдельной работе, более ориентированной на математические аспекты проблемы. Другими словами, ухудшение погрешности сопряжения при кусочной аппроксимации функций с комплексным аргументом должно всегда учитываться, поскольку хорошая точность сопряжения, достигнутая в случае чисто вещественного аргумента, например как в (11), на комплексной плоскости при сопряжении двух кусочно-аппрок-симирующих выражений может стать уже недопустимо большой даже при небольшом удалении от вещественной оси х.
Таким образом, сложное решение ур. (1) о радиальном затухании амплитуды гармонических колебаний в бесконечном напорном пласте с перетеканием на интервале значений Р<1 приближенно можно рассматривать как решение задачи о затухании колебаний в изолированном напорном пласте с несколько измененной эквивалентной пьезопроводностью а*:
£
+ +
о
<3
к
ч
с
§
к
Он
и
0,1
s 0,01
СО
ей
Он
п
0J
W
0,001
Ошибка сопряжения
1,0^ 0,7^ ...
0,5
1,5 2 2,5 3 3,5
Безразмерный аргумент
4,5
5
Р
Рис. 2. Результаты расчета расширенной функции Макдональда с комплексным аргументом в зависимости от безразмерного радиуса: подписи кривых - значения р
(1 + 0,5158Р)2
при 0 < Р< 1.
(18)
На простом физическом уровне полученный результат - выражение (18), легко может быть объяснен. При наличии непроницаемых водоупорных слоев затухание колебаний происходит только плоскорадиально в пределах изолированного пласта. А при перетекании радиальное затухание будет дополнительно усиливаться за счет передачи части энергии колебаний в вертикальном направлении из пласта через слабопроницаемые водоупорные слои. По аналогии с решением задачи о затухании гармонических колебаний в изолированном пласте, такое дополнительное усиление затухания ко-
лебаний в пласте с перетеканием при радиальном удалении от источника колебаний будет эквивалентно уменьшению значения пьезопроводности.
Выводы
Задача о радиальном затухании установившихся гармонических колебаний в напорном пласте с перетеканием может быть эквивалентно аппроксимирована решением аналогичной задачи для случая изолированного пласта с пересчитанным значением эквивалентной пьезопроводности. Пересчитанная эквивалентная пьезопроводность будет меньше фактической, что учитывает дополнительное затухание энергии колебаний при перетекании.
a =
СПИСОК ЛИТЕРАТУРЫ
1. Шестаков В.М. Гидрогеодинамика. - М.: Изд-во МГУ, 1995. -368 с.
2. Синдаловский Л.Н. Справочник аналитических решений для интерпретации результатов опытно-фильтрационных опробований. - СПб.: Изд-во С.-Петерб. ун-та, 2006. - 769 с.
3. Овчинников М.Н. Интерпретация результатов исследований пластов методом фильтрационных волн давления. - Казань: ЗАО «Новое знание», 2003. - 181 с.
4. Бузинов С.Н., Умрихин И.Д. Исследование нефтяных и газовых скважин и пластов. - М.: Недра, 1984. - 269 с.
5. Oldham K.B., Myland J.C., Spanier J. An Atlas of Functions: With Equator, the Atlas Function Calculator. - N.Y.: Springer, 2008. -748 p.
6. Zhang S., Jin J. Computation of special functions. - N.Y.: John Wiley & Sons. Inc, 1996. - 717 p.
7. Люк Ю. Специальные математические функции и их аппроксимации. - М.: Мир, 1980. - 608 с.
8. Математический электронный ресурс для on-line расчета значений специальных функций с высокой точностью. 2011. URL: http://keisan.casio.com/has10/SpecExec.cgi (дата обращения: 27.12.2011).
Поступила 29.03.2012 г.
