Толпаев В.А. Математические модели двумерной фильтрации в анизотропных, неоднородных и многослойных средах

Подождите немного. Документ загружается.

121

а также получаем рекуррентные формулы, выражающие коэффициенты

ksksksks

,,,

−−

βαβα через

1n,k

. Расчётный алгоритм по этим рекуррентным

формулам заключается в следующем.

1. Вычисляются вспомогательные коэффициенты

по формулам:

[]











−−µ++

+−µ+−

=µ

−=

λ−−⋅µ+λ+

λ+−⋅µ+λ−

=µ

−

=µ

−

)k/k1(n/)1n(k/k1

.....................................................................................................

)1n,...,3,2m(;

)1(m/)1m()1(

....................................................................................................

;

1n,k20

n,k

1m,k

m,k

. (12)

В (12) через λ обозначено







=λ

четноеmесли,k/k

нечетноеmесли,k/k

2. Вычисляются вспомогательные коэффициенты

по формулам:

)kk/()kk(

21211k

−=

, а остальные формулы для

получаются из (12) замещени-

ем всех

µ на

(с сохранением указанных в (12) значений индексов).











µ=α

























−













−

µ+













α⋅µ=α













λ−













−

µ+λ+

=α

αµ=α

=α

++−

−+

++−

−+

−

1n,k

kn1n,k

1n,k

1m,k

km1m,k

1m,k

1k2k

q)an(

;

................................................................................

;)am(

;)1(

)1(

................................................................................

;

. (13)

3. Оставшиеся неизвестные

ksksksks

,,,

−−

через заданные коэффициенты

1n,k

и вспомогательные параметры

, находятся из систем уравне-

ний (13) и (14).

122











γ=β

























−













−

γ+













β⋅γ=β













λ−













−

γ+λ+

=β

βγ=β

=β

++−

−+

++−

−+

−

.p)an(

;

................................................................................

;)am(

;)1(

)1(

................................................................................

;

1n,k

kn1n,k

1n,k

1m,k

km1m,k

1m,k

1k2k

. (14)

Порядок решения системы уравнений (13) следующий: из предпоследнего

уравнения по известному

1n,k

находим

, затем по

находим

1n,k −

и т.д. до

α и

α . Затем по значению α

получаем α

−k2

, по

получаем

3k−

α и т.д. до

kn−

Значение же

1n,k +−

находится из последнего уравнения системы (13) через

1n,k

непосредственно. Система (14) решается аналогично.

3.5. Исследование точности фильтрационных расчётов в слоистых средах

методом однородно-анизотропного эквивалентирования

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

слоистая среда заполняет не часть, а всю область фильтрации. Естественно, оцен-

ки точности расчётов фильтрационных потоков, когда анизотропное эквиваленти-

рование применяется ко всей расчётной области, отличаются от тех, когда оно

применяется лишь к части расчётной области. Случай, когда метод однородно-

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

проанализируем на примере расчёта течения жидкости в фильтрационном лотке,

изображённом на рис.19. Для этого снова проводим сопоставительные расчёты

полного фильтрационного потока в слоистой среде на рис.19 и, при тех же гра-

ничных условиях, в анизотропной модели слоистой среды. Применительно к

рис.19 (когда в прямоугольной области фильтрации размещается целое число пар

123

чередующихся изотропных слоёв c одинаковой шириной и с проницаемостями k

и k

) оба метода, локального и интегрального однородно-анизотропного эквива-

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

1мин

kk2

=λ=λ

2макс

=λ=λ , а ГНА всюду направлены по осям x и y декар-

товой системы координат (рис.19) и аналитически задаются функцией (2.5.13) в

виде

(z) = z.

3.5.1. Расчёт полного фильтрационного потока в прямоугольной анизотроп-

ной области

Плоскопараллельное фильтрационное течение жидкости в среде с прямоли-

нейной анизотропией

(z) = z описывается комплексным потенциалом (2.5.14),

который в рассматриваемом случае имеет вид

)(w)y,x(i)y,x(

ϕ , где y/ix

⋅λλ+=ζ . (1)

При этом

⋅λλ

−=ϕ

)y,x(

и )y,x(

, являющиеся соответственно потенциа-

лом скорости фильтрации и функцией тока, применительно к случаю, указанному

на рис.19, должны удовлетворять следующим граничным условиям:

Q;0;

AMDBEC

CDAB

+=ψ=ψ

ρ⋅λλ

−=ϕ=ϕ

, (2)

где Q – полный фильтрационный поток, протекающий в лотке от AB к CD и

подлежащий определению в процессе решения задачи. Из (2) видно, что область

значений w(ζ)- прямоугольник, изображенный на рис.20. Таким образом, для оты-

скания w(ζ), описывающей плоскопараллельную фильтрацию в лотке с пористой

анизотропной средой, необходимо найти функцию, реализующую конформное

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

браженного на рис.19 путем его растяжения вдоль оси y в

/ λλ раз. (Этот полу-

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

плексного переменного

124

Конформное отображение одного прямоугольника на другой осуществлялось

в четыре этапа. На первом этапе нижняя полуплоскость первого вспомогательного

комплексного переменного

(рис.21) отображалась конформно на растянутый

прямоугольник BMDE (рис.19) (в соответствии с аналогичным примером в [78,

123]) при помощи эллиптического интеграла первого рода

)m1)(1(

222

∫

τ−τ−

=ζ

(3)

где

)m(K);m(K/b l= - полный эллиптический интеграл первого рода,

)2/1n(

q21





















∑

∞













⋅

⋅π

−=

expq

. Обращение эллиптического интеграла (3) да-

ет хорошо изученную в математике функцию – эллиптический синус:













ζ⋅

=τ m,

)m(K

. (4)

На втором этапе по формуле (4) вычислялись координаты τ

и τ

точек на

плоскости

, соответствующих значениям

точек A и C на рис.19. На треть-

ем этапе нижняя полуплоскость плоскости

конформно отображалась на нижнюю

полуплоскость второго вспомогательного комплексного переменного t (рис.22)

так, чтобы точки B и C отобразились в точки с координатами t

= -1 и t

=+1, а об-

разы точек A и D отстояли от начала координат на одинаковом расстоянии 1/k, где

0<k<1. Такое отображение осуществляется с помощью дробно-линейной функции,

однозначно определяемой соответствием трех точек A, B и C плоскости

одно-

именным им точкам плоскости t:

k/11

k/1t

+−

−

⋅

−

τ−τ

−

⋅

τ−τ

(5)

125

Неизвестный в последнем уравнении параметр k находится из требования

отображения точки D на рис.21 в одноименную точку на рис.22. Решая получаю-

щееся при этом требовании уравнение, для k найдем

()

−+ρ

, где

τ−τ

−

⋅

τ−τ

−

=ρ

(Для указанного на рис.19 положения

= - 1 и

= 1/m). На четвертом этапе

нижняя полуплоскость плоскости t конформно отображалась на внутренность

прямоугольника на рис.20 при помощи интеграла Шварца-Кристоффеля, имеюще-

го вид (6):

ρ⋅λλ

−

−−

∫

)tk1)(t1(

a)t(w

222

, где

)k(K2

⋅µ⋅

ρ⋅λλ

−=

. (6)

Формула (6) с учетом зависимостей (4) и (5) и определяет искомый ком-

плексный потенциал w(ζ) плоскопараллельного течения в прямоугольной анизо-

тропной области. Для определения с помощью найденного потенциала w(

)

фильтрационного потока Q, заметим, что согласно рис.22 и рис.20

w ⋅=













− .

Вычисляя интеграл (6) при t = - 1/k и сравнивая результат с указанной для w(-

1/k) величиной, найдем поток Q

аниз

(

)

()

kK2

kKgT

аниз

⋅µ

′

⋅ρ⋅λλ

, (7)

где

)k(K

′

- полный эллиптический интеграл первого рода при дополнитель-

ном

модуле

k1k −=

′

. Для частного случая, когда точка A совпадает с M и C сов-

падает с E (рис.19), для величины полного фильтрационного потока потока из (7)

получим значение, равное

l⋅µ

⋅

HgT

. (8)

126

Из (7) и (8) найдём, что безразмерные отношения фильтрационных потоков в

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

(

)

()

аниз

⋅

′

⋅

= . (9)

Формула (9) в вычислительном эксперименте применялась для вычисления

точных значений Q

аниз

для различных отношений

. Результаты расчётов по

формуле (9) отражены в таблице 3.6.

Таблица 3.6

Результаты расчётов полных фильтрационных потоков

мелк

в слоистой среде (рис.19) и в

её анизотропной модели

аниз

. (В вычислительном эксперименте MBED – квадрат; |AB| = |CD|

= 0,5

⋅

H).

Общее число слоев 2 4 10 20 30

Для λ

/ λ

= 0,1 величина Q

/ Q

аниз

= 1,139

/ Q

мелк

1,428 1,338 1,226 1,188 1,167

Погрешность

, %

20,2 17,5 7,1 4,1 2,4

Для

= 0,5 величина Q

/ Q

аниз

= 1,317

/ Q

мелк

1,449 1,397 1,341 1,328 1,321

Погрешность

, %

9,1 5,7 1,8 0,8 0,3

Для

= 0,75 величина Q

/ Q

аниз

= 1,396

/ Q

мелк

1,463 1,433 1,404 1,400 1,396

Погрешность

, %

4,6 2,6 0,6 0,3 0,0

3.5.2. Расчет фильтрационного потока в слоистой области

Аналитический расчёт плоскопараллельного течения в слоистой области на

рис.19 в общем случае встречает серьезные математические трудности. Исключе-

нием служит тот частный случай расположения границ AB и CD, о котором гово-

рилось выше. Если слоистая область содержит в этом частном случае целое число

пар слоёв, то для фильтрационного потока в ней из простых аналитических расчё-

тов получается значение, совпадающее с (8). В общем случае «послойное» реше-

ние задачи сводится к интегрированию системы уравнений Лапласа

127

()

M,....,2,1m,0

y,x

∂

ϕ∂

∂

ϕ∂

=ϕ∆ (10)

для потенциалов

()

⋅

−=ϕ

y,x

всех изотропных слоёв (общее число кото-

рых обозначено через M) с проницаемостями k

. Решения уравнений Лапласа на

границах раздела слоёв Γ

должны удовлетворять условиям сопряжения

1mm

∂

ϕ∂

∂

ϕ∂

, (11)

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

∂

ϕ∂

∂

ϕ∂

∂

ϕ∂

∂

ϕ∂

=−=

, (12)

и на границах CD и AB условиям

⋅

−=ϕ=ϕ

gTk

. (13)

Аналитическое решение задачи (10)-(13) встречает серьёзные математические

трудности. Поэтому эта задача решалась численным методом - методом сеток

[131, 135]. Для этого вся область MBED покрывалась равномерной квадратной

сеткой узлов. Уравнения для потенциалов узлов, расположенных на границах BE,

EC, MD и AM, выводятся из краевых условий на этих участках. Так, например,

для узлов, расположенных на границе BE, выполняется условие

0x/

∂ϕ∂ . Заменяя

производную

x/ ∂ϕ∂ её конечно-разностной аппроксимацией, для )0,i(ϕ получим

,N,...,2,1i);1,i(

)2,i(

)0,i( =ϕ+ϕ−=ϕ (14)

где N - общее число узлов на стороне BE. (Аналогично получаются и записы-

ваются уравнения для узлов, расположенных на границах AM, MD, CE).

Для потенциалов

)j,i(ϕ узлов, лежащих на границах раздела слоёв и узлов,

лежащих внутри изотропных слоёв, пользуемся известными конечно-разностными

представлениями [131, 135]. Значения потенциалов в узлах, расположенных на

границах AB и CD, известны (из граничных условий (13)).

128

Записывая уравнения для потенциалов каждого узла сетки, получим относи-

тельно

)j,i(ϕ систему линейных алгебраических уравнений. Особенность полу-

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

преобразований легко приводится к форме, пригодной для решения методом ите-

раций – применялся метод Зейделя. Достаточные условия сходимости метода Зей-

деля оказываются при этом выполненными. (Выполнение достаточных условий

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

узлов на непроницаемых участках MD, AM, BE, CE подставляли в уравнения для

внутренних узлов области, смежных с граничными). Для уравнений узлов, смеж-

ных с границами AB и CD, условия сходимости метода Зёйделя удовлетворяются

автоматически. Для повышения скорости сходимости итерационного процесса в

вычислениях применяли ускоряющий множитель Либмана [131, 135].

Вычислив значения потенциала в узлах сетки затем можно определить пол-

ный фильтрационный поток Q

мелк

. Для этого в каком-либо m-ом изотропном слое

области фильтрации проводили сечение, параллельное границам MB и DE и вы-

числяли (численным методом по формуле Симпсона) вдоль этого сечения инте-

грал

(

)

()

(

)

мелк

∫∫

−

−+

−

⋅

ϕ−ϕ

=⋅

∂

ϕ∂

, где

(

)

ϕ и

(

)

−

ϕ – потенциалы на линии сече-

ния в выбранном m-ом слое, а h - шаг сетки. Последний интеграл, величина кото-

рого не должна зависеть от выбора места сечения, даёт полный фильтрационный

поток Q

мелк

в слоистой среде. В действительности при расчётах Q

мелк

в качестве

вычисленного значения этой величины бралось среднее арифметическое для 10

равномерно распределённых в области фильтрации сечений. По величине разбро-

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

относительная погрешность расчёта Q

мелк

. Заметим, что во всех расчётах относи-

тельная погрешность Q

мелк

не превышала 2%. Такая точность достигалась за счет

специального выбора шага сетки – проводились последовательные контрольные

расчеты с уменьшением вдвое шага сетки. Расчёты показали, что для достижения

129

в расчётах Q

мелк

указанной точности на стороне BE слоистой области из k слоев

(k=10…30) нужно размещать не менее 5k узлов; из 10 и менее слоев – 50…60 уз-

лов.

Некоторые результаты расчётов потоков Q

мелк

в слоистой области и в её ани-

зотропной модели приведены в таблице 3.6. Общее расположение границ AB и CD

для приведенных в таблице 3.6 расчётов соответствует указанному на рис.19, дли-

ны участков AB и CD брались одинаковыми и равными половине стороны квадра-

та MBED.

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

зотропии, когда

75,0/

=λλ=ε , даже если слоистая среда состоит всего из двух

слоев, она хорошо аппроксимируется анизотропной. При этом относительная по-

грешность δ расчётов потоков

%100

мелк

анизмелк

⋅

−

=δ не превышает 5%. При стрем-

лении

к нулю (когда проницаемости k

и k

изотропных слоёв все резче отлича-

ются друг от друга) погрешность расчётов потоков в слоистых средах, модели-

руемых анизотропными, не будет превышать 5…7%, если: 1) при средней анизо-

тропии, когда

= 0,5 , отношение s/H

≤

0,17 (что соответствует наличию в слоистом

квадрате 3 или более пар изотропных слоев, через s обозначена ширина отдельно-

го слоя); 2) при сильной анизотропии, когда ε = 0,1, отношение s/H≤0,06 (в слой-

чатом квадрате 9 или более пар изотропных полос). При соотношении s/H≤0,03 (в

слойчатом квадрате 15 и более пар слоев) относительная погрешность расчетов

интегральных характеристик методом аппроксимации слоистой среды анизотроп-

ной не превысит 3…5%, если 0,1

≤

Для случаев, когда длины участков AB и CD больше половины стороны

квадрата (большая часть линий тока при этом почти ортогональная границам раз-

дела слоев) точность расчета потоков в слоистой среде методом аппроксимации её

анизотропной оказывается ещё более высокой, чем в таблице 3.6.

130

Если же длины участков AB и CD становятся меньше и меньше, то погреш-

ность расчётов потоков рассматриваемым методом растёт и, начиная с длин около

0,3H, составляет 15…20%.

Выводы по 3-ей главе.

Рассмотренные примеры подтверждают, что со стрем-

лением к нулю отношения a/L (a - толщина отдельного слоя в слоистой среде, L -

характерный размер многослойной области фильтрации) величины таких инте-

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

точностью, аппроксимируя в расчётах слоистые среды их анизотропными моде-

лями. Такая аппроксимация существенно упрощает расчёты и приводит к удобно-

му для анализа аналитическому решению задачи. При этом точность аппроксима-

ции выше у метода интегрального, нежели у метода локального однородно-

анизотропного эквивалентирования. Кроме того, точность метода однородно-

анизотропного эквивалентирования тем выше, чем большая часть линий тока поля

почти ортогональна (или, наоборот, почти параллельна) границам раздела изо-

тропных слоев, составляющих слоистую среду. В проанализированных в 3-ей гла-

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

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

тропной, не превышала 5% при широком диапазоне изменения коэффициента ани-

зотропии (

1/1,0

максмин

<λλ≤

), если a/L

≤

0,03.