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

Выводы ко второй главе.

Во второй главе приведены дополнительные сведения о методе подсеток. Здесь содержится описание построения базовых структур и основных вспомогательных отображений для подсетки, необходимых для эффективного применения метода, поднят вопрос об автоматическом выборе параметров аппроксимации. Доказана теорема о порядке аппроксимации потока при использовании метода подсеток в изотропной однородной среде. от переменных давления в блоках шаблона и иметь ненулевые производные по ним, которые входят в строку матрицы якобиана в методе Ньютона (см. раздел 1.1). При использовании многоточечного метода размер этой матрицы не меняется, но растет ширина ее шаблона1 в соответствии с шаблоном метода, что напрямую влияет на скорость сходимости неявной схемы. Как правило для оценки величины шаблона метода используют количество элементов в нем для внутреннего блока квазиравномерной сетки. Например, О-метод имеет 27-точечный шаблон, задействуя тем самым все элементы в окрестности блока, а производный от метода опорных операторов метод базисов связей — 19-точечный шаблон, [7]. Для метода подсеток шаблон состоит из максимального количества — 27 точек.

Для многих многоточечных методов происходит так называемое вырождение шаблона в равномерном случае при диагональном тензоре проницаемости. Это означает, что значения давления в некоторых блоках, которые входили в шаблон, перестают влиять на давление в рассматриваемом блоке. Обычно шаблон вырождается в 7-точечный шаблон ТРБА. Это происходит для Ь-метода, О-метода, метода базисов связей и других. С одной стороны это полезное свойство, поскольку оно позволяет уменьшить количество ненулевых коэффициентов для тех строк якобиана, которые соответствуют вырожденному шаблону, а значит, и увеличить скорость расчета. С другой стороны это означает, что метод сильнее подвержен ориентационному эффекту (см. обзор литературы). Другими словами, отсутствует течение «по диагонали» в случае, близком к равномерному. Если даже окажется, что давление в одном из блоков, лежащих по диагонали, сильно упадет по сравнению с давлением в нашем блоке, то при вырождении шаблона потоки через интерфейсы центрального блока останутся прежними. При больших по величине временных шагах схемы данное обстоятельство может вносить значимую ошибку (например, в некоторых случаях, когда в набор данных перестают добавлять новые.

В3 Е.

В1.

В 2.

Рис. 3.15. Случай нефизичного течения для метода базисов связей перфорации, временной шаг может достигать нескольких месяцев).

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

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

Пусть тензор проницаемости единичный и давление в блоке В{ равно Рг. Приближение градиента в центре Вг будет равно ¦ ' где — расстояние между центрами блоков Вь и В^. Нормаль к грани Е, внешняя по отношению к В, равна, а > 0. Это означает, что приближение потока, построенное по этому базису равно.

Если окажется, что рз > р и р3 — р < — р2), то значение потока.

О 3 о 2.

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

Отметим еще одну особенность метода подсеток, касающуюся взаимных перетоков между соседними блоками. Рассмотрим набор данных с сеткой из четырех блоков Ь, Ь2, 64 размером 2x2x2 (см. рис. 3.17), для каждого из них зададим тензор проницаемости К, К2, К2, К соответственно, где.

Л л Л.

1 (Л.

Кх V.

1 з о 0 0 1.

Ко =.

Л 0 0.

0 Л 0.

0 0 1.

На границе резервуара поставим условия непротекания. Вычислим потоки при применении О-метода аналитически. Резервуар распадается на девять областей взаимодействия: четыре угловых, которые не используются для вычисления потоков, четыре приграничных, которые похожи в силу симметрии, и.

Рис. 3.18. Области взаимодействия для тестового набора данных.

VI 41 Р2.

42 7 777 777 777 777 778 282 705 982 024 253 779 738 624 д’л Рис. 3.19. Область взаимодействия первого типа для тестовой набора данных одну центральную (см. рис. 3.18). Итак, нам нужно рассмотреть всего два типа областей.

В первом случае в область входят два блока Ъ2. Обозначим давления в их центрах через р и р2. При использовании О-метода предполагается, что давление в каждом блоке линейно, непрерывно в центрах интерфейсов, и имеет место равенство потоков через части интерфейсов, попавших в область взаимодействия. Обозначим давления в центрах интерфейсов через д2, <7з (см. рис. 3.19). Тогда бі = (01 ~РЪ Р1- (/2)Т, Кі^7и~Ьі = (3.

С.

— с.

V1 91 Р2.

Рис. 3.20. Область взаимодействия второго типа для тестового набора данных.

Равенство потоков имеет вид.

3¹ — Ф2 — 2р1 = Л (р2 — 91), — Зд2 + 2рх = 0, ЧР2 — Яз) = о.

Решив систему, получим давление.

8¹ + 3 Хр2 91 ~ 8 + ЗА и поток через половину интерфейса между Ъ и 1)2, равный.

8Л / 3 = НР2 — Я) = (Р2 — Р1)8 зл = (Р2 — Р1) Л (1 + -л (Р2 -^1)Л + 0(Л2), при, А —>> 0.

Во втором случае в область входят все четыре блока из набора данных. Обозначим давления в них через р, р2, рз, Р4, а давления в центрах интерфейсов через Яз, Яа (см. рис. 3.20). Рассмотрим только те уравнения, которые необходимы для вычисления потока через интерфейс между блоками свойствами при малом Л потоки вычисляются по следующим формулам:

2 = —4.312 839/?! + 2.830 662² — 1.290 579р3 + 2.772 756р4, /з = -4.312 839рі - 1.290 579² + 2.830 662р3 + 2.772 756/?4, /24 = —2.772 756рі - 2.830 662² + 1.290 579р3 + 4.312 839/?4, /з4 = —2.772 756рі + 1.290 579р2 — 2.830 662/?3 + 4.312 839р4, где — поток через интерфейс между блоками Ьі и в направлении последнего. Заметим, что при // > 0 жидкость перетекает из Ъ^ в Ьі, поскольку ее движение происходит в противоположном направлении относительно градиента (от большего потенциала к меньшему). Общий переток через все интерфейсы по направлению от блока &і к блоку 64 равен = -(/і + Л3 + /2 + /з) = 14.17 1190(рір4). В то же время отток жидкости из блока Ь2 равен -1.540 083/?! + 5.661 324/?2 — 2.581 158/?3 — 1.540 083/э4.

В начальный момент времени при давлениях (3.1) имеем равенства, а = 708.5595, Ъ = 77.0042, т. е. скорость течения «по диагонали» на порядок выше, чем изменения в других блоках со слабой проводимостью. Для ¿—метода на этом же наборе данных получились значения а= 168.2964, Ь — 84.1482.

Полученный результат подтверждается на практике. На рис. 3.21 приведены результаты расчетов для описанного набора данных из четырех блоков.

50.000.

100.000.

Рис. 3.21. Карты давления (атм.) для тестового набора данных при расчете различными методамипо центру — начальное давление, слева — конечное давление при расчете ¿—методом, справа — конечное давление при расчете методом подсеток.

Видно, что для метода подсеток диагональный переток значительно больше, чем при применении Ь-метода.

Попытаемся разобраться, является ли такое поведение модели при расчете методом подсеток адекватным, построив специальный тестовый набор данных. При повороте на 45 градусов резервуара, изображенного на рис. 3.17, тензор проницаемости во всех блоках становится диагональным (К2 является скалярным тензором в Оху, а оси К[ совпали с осями новых координат). Для того, чтобы мы могли рассчитать набор данных двухточечной схемой, наложим на резервуар равномерную ортогональную сетку с малым шагом, получив тем самым разбиение крупных блоков на мелкие. Будем обозначать области нового набора данных, соответствующие этим блокам изначальных данных, через Ь{ при 1 < г ^ 4. Зададим начальное давление (см. рис. 3.22) и тензор проницаемости (см. рис. 3.23) схожим с предыдущим примером способом. Матрица К2 осталась прежней при повороте, а стала равной.

В центральном блоке величина тензора была специально уменьшена, чтобы сократить прямую связь областей ?1 и 64 с высокой проводимостью. На.

4 0 0 Л.

К[ = 0 2 0 альным искусственным материалом (называемым проппантом), разрывающая нефтяные пласты резервуара. С математической точки зрения эта процедура локально изменяет некоторые параметры набора данных, в том числе тензор абсолютной проницаемости. Может оказаться, что полученная с помощью ГРП трещина не проходит вдоль координатных осей, а имеет сложную структуру. Проводимость вдоль трещины бывает в 30−50 раз выше, чем проводимость самого пласта. Это вызывает сильную анизотропию с недиагональным тензором проницаемости. Продемонстрируем это на примере.

Рассмотрим два преобразования тензора проницаемости, которые нам потребуются для описания трещин. Введем матрицу поворота на угол <�р> относительно оси Ог cos <�р — sin ср (Л.

А = А{<�р) sin (f cos ip О 0 0 1 Если М = AKAr, то выполнены равенства V.

Ми.

М22 М33 Ми Mi3 М23.

Ки COS2 ср — 2К2 sin (р cos ср + К22 sin2 Кц sin2 (р + 2К2 sin COS ip + К22 cos2 (р, Кзз:

Кц — К22) sin ip cos (р + К12(cos2 (р — sin2 ср), Ki3 COS ср — К23 sin (р, Кsin (р + К2з cos (р.

Введем матрицу поворота на угол ф относительно оси Оу:

В = В (яр) = cos ф 0 — sin ф^ 0 10^sin ф 0 cos ф у.

Если М = ВКВТ, то выполнены равенства.

Мц = #11 COS2 ф — 2#i3 sin ф cos ф + #33 sin2 ф,.

М22 = #22,.

Мзз = #11 sin2 ф + 2#i3 sin ф cos ф + #33 cos2 ф, М12 = #12 COS Ф — #23 sin Ф, m13 = (#11 — #33) sin cos Ф + #13 (cos2 ф — sin2 ф), M23 = #12 sin Ф + #23 cos ф.

Трещина в направлении Ох задается таким диагональным тензором п О О.

О #22 о V о.

О # зз у для которого #11 #22 и #п #33. Тензор К' = при = arctg будет описывать трещину, идущую к блоку, лежащему по диагонали в том же слое равномерной сетки (здесь х, у, г — размеры блока-параллелепипеда). Учитывая диагональность #, имеем j и —.

К'22 =.

К*.

К[2 = о о.

11 COS (р + #22 sin </?,.

11 sin2 У? + #22 COS2 (р, (#11 — #22) sin у? COS <�р,.

О, 0.

23 =.

Тензор К" = АВКВТАт при у? = arctg ^ и ф = аг^, х будем использовать для трещины, идущей к блоку, лежащему по направлению главной диагонали в соседнем слое равномерной сетки. Учитывая диагональность #,.

1 ¦¦

Г с.

Мг я, а я мая яял.

I шшш тятя.

Я ятта.

1 II 1 I о. оо мнннмна 1 .оо.

Рис. 3.28. Карта насыщенности нефти во втором слое в конце 40-го временного шага на тестовом наборе данных без трещин при расчете двухточечным методом.

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

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

28. Aavatsmark I., Barkve T., Boe O., Mannseth T. Discretization on Unstructured Grids For Inhomogeneous, Anisotropic Media. Part II: Discussion And Numerical Results // SIAM Journal on Scientific Computing. 1998. Vol. 19, no. 5. Pp. 1717−1736.

29. Aavatsmark I., Eigestad G., Heimsund B.-O. et al. A New Finite-Volume Approach to Efficient Discretization on Challenging Grids // SPE Journal. 2010. Vol. 15, no. 3. Pp. 658−669.

30. Aavatsmark I., Eigestad G., Mallison B., Nordbotten J. A compact multipoint flux approximation method with improved robustness // Numerical Methods for Partial Differential Equations. 2008. Vol. 24, no. 5. Pp. 1329−1360.

31. Andreev A. B., Dautov R. Z. Boundary-flux estimates for a finite element approximations // Godishnik Vissh. Uchebn. Zaved. Prilozhna Mat. 1989. Vol. 25, no. 4. Pp. 105−113.

32. Au A. D. K., Behie G. A., Rubin B., Vinsome P. K. Techniques for Fully Implicit Reservoir Simulation // SPE paper 9302. Presented at SPE Annual Technical Conference and Exhibition, Dallas, Texas. 1980. 12 pp.

33. Aziz K., Settari A. Petroleum reservoir simulation. London: Applied Science Publishers, 1979.

34. Brand C., Heinemann J., Aziz K. The Grid Orientation Effect in Reservoir Simulation // SPE paper 21 228. Presented at SPE Symposium on Reservoir Simulation, Anaheim, California. 1991. 12 pp.

35. Brezzi F., Fortin M. Mixed and hybrid finite element methods. Springer series in computational mathematics. New York: Springer-Verlag, 1991.

36. Chen Z. Reservoir Simulation: Mathematical Techniques in Oil Recovery. Society for Industrial and Applied Mathematics, Philadelphia: SLAM, 2007.

37. Chen Z., Huan G., Ma Y. Computational Methods for Multiphase Flows in Porous Media. Philadephia, PA: SIAM, 2006.

38. Coats K. H. IMPES Stability: Selection of Stable Timesteps // SPE Journal. 2003. Vol. 8, no. 2. Pp. 181−187.

39. Craft B., Hawkins M. F. Applied Petroleum Reservoir Engineering. 2 edition. NJ 7 632: Prentice Hall PTR Englewood Cliffs, 1991.

40. Danilov A. A., Vassilevski Y. V. A monotone nonlinear finite volume method for diffusion equations on conformal polyhedral meshes // Russ. J. Numer. Anal. Math. Modelling. 2009. Vol. 24. Pp. 207−227.

41. Economides M. J., Hill A. D., Ehlig-Economides C. Petroleum Production Systems. Upper Saddle River, New Jersey: Prentice Hall Ptr, 1993.

42. Ertekin T., Abou-Kassem J. H., King G. R. Basic Applied Reservoir Simulation. Richardson, Texas: Society of Petroleum Engineers, 2001.

43. Fanchi J. R. Principles of Applied Reservoir Simulation. 3 edition. Atlanta, GA: Elsevier Science Ltd, 2005.

44. Gunasekera D., Childs P., Herring J., Cox J. A Multi-Point Flux Discretization Scheme for General Polyhedral Grids // SPE paper 48 855. Presented at SPE International Oil and Gas Conference and Exhibition, Beijing, China. 1998. 9 pp.

45. Hoteit H., Mose R., Philippe B. et al. The maximum principle violations of the mixed-hybrid finite-element method applied to diffusion equations // Interna.

54. Verma S., Aziz K. A Control Volume Scheme for Flexible Grids in Reservoir Simulation // SPE paper 37 999. Presented at SPE Reservoir Simulation Symposium, Dallas, Texas. 1997. 13 pp.

55. Walsh M. P., Lake L. W. A Generalized Approach to Primary Hydrocarbon Recovery (Handbook of Petroleum Exploration and Production). 4 edition. Atlanta, GA: Elsevier Science Ltd, 2003.