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

Механико-математический факультет Московский государственный университет им. М.В. Ломоносова

Мельниченко Никита Сергеевич

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

05.13.18 - Математическое моделирование, численные методы и комплексы

программ

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

( 1 СЕН 2011

Москва-2011

4852509

Работа выполнена на механико-математическом факультете Московского государственного университета им. М. В. Ломоносова.

Защита состоится «22» сентября 2011 г. в 16 часов 00 минут на заседании диссертационного совета Д 212.081.21 при ФГАОУВПО «Казанский (Приволжский) федеральный университет» по адресу: 420008, г. Казань, ул. Кремлевская, д. 18, корп. 2, ауд. 218.

С диссертацией можно ознакомиться в научной библиотеке им. Н. И. Лобачевского ФГАОУВПО «Казанский (Приволжский) федеральный университет».

Ведущая организация:

Официальные оппоненты:

Научный руководитель:

кандидат физико-математических наук, доцент,

Богачев Кирилл Юрьевич

доктор физико-математических наук,

доцент,

Федотов Евгений Михайлович кандидат физико-математических наук, Капырин Иван Викторович Московский энергетический институт (технический университет)

Автореферат разослан «.

№ » о^оЩ^л 2011 г.

Ученый секретарь

диссертационного совета,

доктор физико-математических наук

Задворнов О. А.

Общая характеристика работы

Актуальность работы. С ростом производительности современных вычислительных систем увеличивается сложность гидродинамических моделей, используемых при анализе месторождений. В настоящее время начинает пре-^бладатьточкГзрШйя7чтб"расчё^ в резервуаре надо произ-

водить на геологических моделях, чтобы избежать потерь данных при укрупнении сетки (ирБсаНвд). Описание сложных геологических структур (таких как разломы, выклинивание и др.) может потребовать использования неструктурированных сеток.

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

Для сокращения погрешности используют различные многоточечные методы (МРРА) — базисов связей, О-метод, II-метод и др. — которые для вычисления потока через грань между двумя блоками сетки привлекают значения в соседних к ним блоках. Эти методы хорошо приближают поток, не ортогональный грани блока, но не решают проблему ориентационного эффекта и добавляют еще одну проблему — возникновение нефизичных потоков из блока с меньшим давлением в блок с большим давлением при сильной неортогональности граней координатным осям и сильной анизотропии тензора абсолютной проницаемости.

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

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

Цели диссертационной работы:

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

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

3. Проверить адекватность метода на тестовых и реальных наборах данных геологического масштаба, а также провести его сравнение с другими методами аппроксимации.

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

• при использовании метода подсеток нефизичные перетоки практически отсутствуют,

• метод подсеток снижает ориентационный эффект схемы.

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

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

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

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

На защиту выносятся следующие основные результаты и положения:

1. Предложен новый метод многоточечной аппроксимации потока для задач фильтрации — метод подсеток.

2. Приведено теоретическое обоснование качества аппроксимации потока для изотропной однородной среды.

3. Разработан программный модуль, реализующий данный метод, для промышленного комплекса гидродинамического моделирования, и произведена проверка адекватности метода на основе его сравнительного анализа с другими методами МРРА и согласованности численных экспериментов с физической картиной явления.

Апробация работы. Основные результаты диссертации докладывались на различных конференциях и семинарах. Среди них 4-я международная конференция «Математические идеи П.Л. Чебышева и их приложения к современным проблемам естествознания» (Обнинск, 2008 г.), конференция «Ломоносовские чтения» (Москва, 2008, 2009 гг.), научный семинар кафедры вычислительной математики механико-математического факультета МГУ им. М.В. Ломоносова (Москва, 2008 г.), научный семинар Института математического моделирования РАН (Москва, 2008 г.), научный семинар кафедры дифференциальных уравнений Московского энергетического института (Москва, 2009 г.), научный семинар кафедры прикладной математики МГТУ им Н.Э. Баумана (Москва, 2010 г.), научный семинар кафедры вычислительной математики Казанского федерального университета (Казань, 2011 г.)

Публикации. По теме диссертации опубликовано 5 статей в реферируемых журналах [1-5], четыре из которых рекомендованы ВАК РФ для публикации результатов кандидатских диссертаций.

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

Структура и объем диссертации. Диссертация состоит из введения, обзора литературы, трех глав, заключения и библиографии. Общий объем дис-

сертации 124 страницы, из них 117 страниц текста, включая 52 рисунка и 1 таблицу. Библиография включает 55 наименований на 7 страницах.

Содержание работы

Во Введении обоснована актуальность диссертационной работы, сформулирована проблема, описывается структура диссертации.

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

где пс — количество компонентов смеси; i?gj0 — растворимость газа в нефтяной фазе; Bi — коэффициент объемного расширения фазы 1,1 — о (oil), g (gas), w (water); ф — пористость среды; xc,i — молярная доля компонента с в фазе 1\ — молярная плотность фазы I; к — тензор абсолютной проницаемости; кг[ — относительная проницаемость фазы I; щ — вязкость фазы /; 7; — величина вертикального градиента гидростатического давления в фазе l\d — функция глубины; Рсоg — капиллярное давление в системе нефть-газ; Pcow — капиллярное давление в системе вода-нефть; qc — источник компонента с (скважина).

Неизвестными в этой системе являются:

1. Nc = Nc(t, х, у, z) — молярная плотность компонента с (для модели черной нефти компонентами служат вода, нефть и газ);

jt {фЫс) = div Y, ~ 7/Vd)) +

д_ dt

+ Чс

(1)

2. 5; = ж, у, г) — насыщенность фазы I;

3. р1 = р1(£, х, у, г) — давление в фазе I.

Для системы уравнений (1) задаются начальные условия, а также на внешней границе резервуара ставятся однородные условия Неймана (непротекания).

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

Для иллюстрации нефизичного

течения рассмотрим систему из трех блоков (см. рис. 1). Пусть тензор проницаемости единичный и давление в блоке В1 равно р,. Нормаль к грани Е, внешняя по отношению к В\, рав-

на

(«л)

Л*2+1

, а > 0. Приближение потока

в центре Въ построенное по методу

базисов связей равно

1 / Р2 - Р\ , РЗ ~ Р\ ос-:--Ь

Рис. 1. Случай нефизичного течения для мсто-г/^ПГ ¿12 ' <*13 У' Да базисов связей где ¿у — расстояние между центрами блоков В; и В]. Если окажется, что р.з > р\ и рз — Р1 < '-^-{рг — Р2), то значение потока станет отрицательным, т.е. будет переток из В\ в В3 (хотя рз > Рг > Р2). увеличивающийся тем сильнее, чем меньше давление Р2- В анизотропной среде для блоков, размеры которых отличаются в десятки раз по различным направлениям, этот эффект проявляется еще сильнее.

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

Нефизичные перетоки присутствуют для других многоточечных методов. Каждый из методов имеет область применимости, вне которой использование метода нежелательно, поскольку это приводит к нефизичным эффектам. Среди остальных методов ¿-метод имеет наибольшую область применимости, однако на реальных грубых сетках он довольно часто приводит к аномальному течению. Таким образом, появление нефизичных перетоков является одной из главных проблем многоточечных методов.

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

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

На плоских ортогональных сетках ориентационный эффект как правило выражается в отсутствии течения «по диагонали». На рис. 2 изображены блоки, значения в которых влияют на дискретный поток через интерфейсы центрального блока, в случаях с ортогональной сеткой при разбиении вдоль координатных осей (сплошной линией) и с такой же сеткой, повернутой на 45° (пунктиром). Рис. 2. Схема перетоков при наличии ориента- „

^ Видно, что в первом случае переток

ционного эффекта

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

В первой главе подробно изложен многоточечный метод подсеток.

В разделе 1.1 приведена пространственная аппроксимация уравнений и отмечен общий ход решения задачи. Сетка представляется в виде набора множеств узлов ЛГ, блоков В и интерфейсов Ь, через которые происходят перетоки. Введены понятия шаблона метода относительно интерфейса и коэффициентов проводимости t^:, которые определяют действие дискретного аналога

оператора J(/CV-,n) ds на сеточную функцию (здесь I е L, а матрица К со-i

держит тензор абсолютной проницаемости к). В случае, когда в шаблон входит ровно две точки для всех интерфейсов I, метод аппроксимации называется двухточечным (TPFA, Two Point Flux Approximation), а когда существуют интерфейсы с более чем двумя точками в шаблоне — многоточечным методом (MPFA, Multipoint Flux Approximation).

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

В разделе 1.2 описывается построение подсетей внутри так называемой области взаимодействия вокруг узла сетки. Она представима в виде набора базовых тетраэдров PMEV (см. рис. 3), где Е — центр содержащего узел V ребра, М — центр содержащего его интерфейса, Р — центр содержащего его блока.

Построение подсетей в области взаимодействия осуществляется пу- Рис 3 область взаимодействия в трсхмсриом тем разбиения базовых тетраэдров на случае более мелкие, которые будут являться

элементами подсетки. Для этого задается единый параметр разбиения т] > 2, который означает количество узлов подсетки на ребре базового тетраэдра. Вначале строится множество узлов разбиения тетраэдра (см. рис. 4, а). Их

3 2

количество составляет N(7]) = ^ + % + Затем базовый тетраэдр разбивается на маленькие элементы-тетраэдры трех типов. Элементы первого типа (см. рис. 4, б) подобны базовому и их количество равно N(7] — 1). Второй тип элементов (см. рис. 4, е) появляется уже при ту = 3. В этом случае при отбрасывании всех четырех элементов первого типа (углы тетраэдра) остается октаэдр, который разделяется внутренним отрезком на 4 тетраэдра. При увеличении г] множество таких октаэдров растет, а общее количество полученных из них элементов второго типа равно АЫ{г] — 2). Элементы третьего типа (см. рис. 4, г) появляются при г/ ^ 4. Они находятся в окружении октаэдров, являются подобными базовому, их количество равно N(17 — 3). Таким образом, число маленьких тетраэдров при разбиении одного базового составляет

(V -1)\

Описанное разбиение устроено таким образом, что все грани базовых тетраэдров разбиты одинаково на (77 — I)2 треугольников. Это позволяет рассматривать подсетку как структурированную сетку из тетраэдров на всей области взаимодействия.

В разделе 1.3 представлены различные варианты краевых задач на под-сетке. Рассматривается задача о нахождении приближения ик функции давления в конечномерном пространстве и11 непрерывных функций, линейных в каждом элементе подсетки. Это позволит в дальнейшем вычислить поток через части интерфейсов сетки, попавших в область взаимодействия. Для этого сначала находятся приближенные значения в граничных узлах подсетки, а затем вычисляется дискретный аналог решения задачи сНу КУи = 0.

Особое внимание в этом разделе уделено составлению граничных условий задачи. Граница области взаимодействия может быть представлена в виде

Рис. 4. Разбиение базового тетраэдра при т] = 4

объединения поверхности Гь принадлежащей границе резервуара, и поверхности Гг, проходящей внутри расчетной области. На 1\ выбираются граничные условия резервуара, а на Гг специальным образом ставятся условия первого рода, основываясь на значениях функции в центрах блоков и граничных условиях на соседних интерфейсах сетки.

Обозначим рх,р3 давления в блоках сетки, попавших в область взаимодействия, величины давления на интерфейсах первого рода через -ф\,..., фб1, а величины потока на интерфейсах второго рода — через ф\,...,ф82 (значения 5, 5ь .52 зависят от области взаимодействия). Тогда значения в граничном узле

подсетей А представляются в виде

и\

(Л) = а04Рг + ^ 01,¿^г + ^ 0,2,гф,

(2)

¿=1

{=1

где ао,и Я1,ь «2,; — коэффициенты, зависящие только от геометрии области и тензора проницаемости К.

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

Для его получения необходимо решить (в + вх + ¿г) систем линейных уравнений с одинаковой положительно определенной симметричной матрицей А. В разделе 1.3.2 детально описан метод вычисления А и правых частей систем по данным на сетке и подсетке. Сами системы предлагается решать методом сопряженных градиентов с предобуславливателем неполная факторизация Хо-лецкого.

В разделе 1.4 приведен алгоритм вычисления коэффициентов проводимости и шаблона, основываясь на решениях (3). Отмечено, что составление и решение задач на подсетке является независимым для различных узлов сетки, и эффективность распараллеливания метода подсеток будет высокой.

Во второй главе изложены способы вычисления параметров метода подсеток, а также доказана теорема о порядке аппроксимации потока. Проанализированы различные аспекты реализации метода, включая вычисления на параллельной ЭВМ.

В разделе 2.1 подробно изложено построение различных структур и отображений для подсетей, которое было опущено в предыдущей главе. Описывается алгоритм создания подсетей в такой же форме, как и сетка — в виде

(3)

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

-В-разделе-272-описан-один-из-способов-приближения-градиента-функ~

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

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

В разделе 2.4 доказана теорема о порядке аппроксимации потока в однородном изотропном случае. Предположим, что в некоторый момент времени мы точно знаем значения функции и в центрах блоков, а диаметр области взаимодействия Q равен h. Будем говорить, что векторы v¡,..., vn удовлетворяют Я-критерию, если

п

(vi,...,vn) > оПЫ, ¿=1

Ы > С2/1, 1 ^ i ^ п,

где (vi,...,v„) — объем n-мерного параллелепипеда, натянутого на векторы, ci, С2 > 0 — фиксированные константы. При этом тетраэдр ABCD удовлетворяет Я-критерию, если векторы AÉ, AÓ, Añ удовлетворяют Я-критерию.

Теорема 1. Пусть все базовые тетраэдры области взаимодействия Í! удовлетворяют Я-критерию. Тогда порядок аппроксимации потока через часть а интерфейса, попавшую в область взаимодействия внутреннего узла сетки, при расчете с помощью метода подсеток равен mes(er)0(/i) для функции и е С2{П).

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

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

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

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

(а). Набор данных с анизотропным тензо- (б). Срезы множества блоков с насыщенно-ром, повернутым на 45° стью нефти менее 0.9

0.00 • ■■■■■ 1.00 Рис. 6. Недиагональный тензор

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

(а). Без трещин (б). С диагональной трещиной

0.00 ^ЗИНН 1.00

Рис. 8. Нефтенасыщенность для тестового набора данных с ГРП

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

В разделе 3.2 рассмотрены тесты с недиагональным тензором проницаемости. Для набора данных с ортогональной сеткой и тензором проницаемости, который^вМется^иагоналБным-в-гистеме-координа^-повернутой-на-45°, фронт принял форму эллипса, повернутого на тот же угол (см. рис. 6, а). Исследования показали, что использование сеток различной сложности (с локальным сгущением, сдвигами и пр.) вызывают незначительные оптические искажения фронта и также дают качественно правильный результат (иллюстрации приведены в работе). Для полноразмерного набора данных в однородной анизотропной среде фронт имеет форму эллипсоида (см. рис. 6, б), что согласуется с теорией. Таким образом, метод подееток определяет качественно правильное поведение модели.

В разделе 3.3 приведено сравнение метода подееток с другими многоточечными методами. Как правило для оценки сложности метода используют количество элементов в шаблоне относительно внутреннего блока квазиравномерной сетки (иначе говоря, минимальному количеству значений давления в некоторых блоках сетки, достаточному для определения потока через все грани рассматриваемого блока). Для метода подееток шаблон состоит, как и для О-метода, из максимального количества — 27 точек. Кроме того, этот шаблон не подвержен вырождению на ортогональной сетке с диагональным тензором, в отличии от шаблонов остальных рассмотренных методов. Автоматическое вырождение позволяет повысить скорость расчета в некоторых случаях, однако уменьшение шаблона метода означает, что метод плохо справляется с ориентационным эффектом, поскольку отсутствует течение «по диагонали» в равномерном случае. Метод подееток всегда имеет полный шаблон, но за счет этого позволяет снизить ориентационный эффект, вносимый комбинацией используемых сетки и тензора. Это хорошо видно на следующем примере.

Рассмотрим специальный набор данных с сеткой из четырех блоков с таким начальным давлением, как показано на рис. 7. Для левого нижнего и правого верхнего блоков зададим недиагональный тензор, направленный по диагонали, а для двух других — оставим скалярный, но уменьшим его по величине на два порядка. В этом случае расчет О-методом получается по характеру очень похожим на расчет двухточечной схемой — переток идет только через блоки с низкой проводимостью. В то же время при расчете методом подсеток появляется течение «по диагонали». Чтобы проверить адекватность такого поведения, повернем систему координат на 45°. Матрица тензора в каждом блоке при этом становится диагональной, и мы можем выполнять расчет двухточечной схемой, если наложим на резервуар мелкую ортогональную сетку. Полученные в работе результаты для такого набора данных говорят в пользу метода подсеток, показывая, что наблюдается течение между блоками с высокой проводимостью.

Описанное свойство метода подсеток снижать ориентационный эффект нашло свое применение при моделировании трещин, возникающих при гидроразрыве пласта (ГРП), который широко применяется в настоящее время для увеличения нефтеотдачи действующих добывающих скважин. Может оказаться, что полученная с помощью ГРП трещина не проходит вдоль координатных осей, а имеет сложную структуру. Проводимость вдоль трещины бывает в 30-50 раз выше, чем проводимость самого пласта. Это вызывает сильную анизотропию с недиагональным тензором проницаемости. На рис. 8 изображены карты насыщенностей для набора данных с диагональной трещиной и без нее. Видно, что прорыв воды при наличии ГРП происходит раньше, что может существенно сказаться на режиме работы скважины.

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

абсолютной проницаемости показало, что метод подсеток значительно снижает количество групп блоков, в которых появляется нефизичные перетоки. Его использование совместно с ¿-методом и упрощенным методом подсеток, когда для задач на подсетке ставятся условия непротекания, по специальной методике, описанной в раЬоте, позволяётпрактическитобавитьсятгганомальных-течений.

В Заключении сформулированы основные результаты диссертации.

Основные результаты работы

1. Предложен новый метод многоточечной аппроксимации потока для задач фильтрации — метод подсеток.

2. Приведено теоретическое обоснование качества аппроксимации потока для изотропной однородной среды.

3. Разработан программный модуль, реализующий данный метод, для промышленного комплекса гидродинамического моделирования, и произведена проверка адекватности метода на основе его сравнительного анализа с другими методами МРБА и согласованности численных экспериментов с физической картиной явления.

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

Список публикаций

1. Богачев К. Ю., Мельниченко Н. С., Шелков В. Г. Применение метода опорных операторов для эффективного гидродинамического моделирования пластов на неструктурированных неортогональных сетках // Нефтяное хозяйство. 2008. № 2. С. 58-59.

2. Богачев К. Ю., Мельниченко Н. С. О пространственной аппроксимации методом подсеток для задачи фильтрации вязкой сжимаемой жидкости в пористой среде // Вычислительные методы и программирование. 2008. Т. 9. С. 191-199.

3. Богачев К. Ю., Мельниченко Н. С., Шелков В. Г. Применение метода подсеток для гидродинамического моделирования в анизотропной среде на нерегулярных сетках // Нефтех. 2008. № 4. С. 8-12.

4. Мельниченко Н. С. Граничные условия локальной задачи в методе подсеток // Программные продукты и системы. 2009. № 3. С. 11-13.

5. Мельниченко Н. С. Об адаптивной аппроксимации задачи фильтрации вязкой сжимаемой жидкости в трещиноватой пористой среде // Вестник Московского Университета. Серия 1: Математика, Механика. 2009. № 6. С. 51-52.

Заказ № 70-А/08/2011 Подписано в печать 03.08.2011 Тираж 100 экз. Усл. п.л. 1

ООО "Цифровичок", тел. (495) 649-83-30 > \vw\v. с/г. ги; е-таН: т/о@с/г. ги

Введение.

Обзор литературы.

Глава 1. Аппроксимация методом подсеток.

1.1. Многоточечная аппроксимация задачи фильтрации.

1.2. Построение подсетей.

1.3. Краевые задачи на подсетке.

1.4. Вычисление коэффициентов проводимости.

Глава 2. Вычисление параметров метода подсеток.

2.1. Построение структур и отображений для подсетей.

2.2. Задача приближения градиента функции.

2.3. Вычисление значений функции на границе подсетей.

2.4. Порядок аппроксимации потока.

Глава 3. Численные результаты.

3.1. Сравнение с двухточечной аппроксимацией.

3.2. Использование недиагонального тензора проницаемости

3.3. Сравнение с другими методами многоточечной аппроксимации

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

Во второй главе приведены дополнительные сведения о методе подсеток. Здесь содержится описание построения базовых структур и основных вспомогательных отображений для подсетки, необходимых для эффективного применения метода, поднят вопрос об автоматическом выборе параметров аппроксимации. Доказана теорема о порядке аппроксимации потока при использовании метода подсеток в изотропной однородной среде. от переменных давления в блоках шаблона и иметь ненулевые производные по ним, которые входят в строку матрицы якобиана в методе Ньютона (см. раздел 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 7777777777777777777777777777777777777777 д'л Рис. 3.19. Область взаимодействия первого типа для тестовой набора данных одну центральную (см. рис. 3.18). Итак, нам нужно рассмотреть всего два типа областей.

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

С

- с

V1 91 Р2

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

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

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

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

8^1 + 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.312839/?! + 2.830662^2 - 1.290579р3 + 2.772756р4, /з = -4.312839рі - 1.290579^2 + 2.830662р3 + 2.772756/?4, /24 = —2.772756рі - 2.830662^2 + 1.290579р3 + 4.312839/?4, /з4 = —2.772756рі + 1.290579р2 - 2.830662/?3 + 4.312839р4, где — поток через интерфейс между блоками Ьі и в направлении последнего. Заметим, что при // > 0 жидкость перетекает из Ъ^ в Ьі, поскольку ее движение происходит в противоположном направлении относительно градиента (от большего потенциала к меньшему). Общий переток через все интерфейсы по направлению от блока &і к блоку 64 равен = -(/і + Л3 + /2 + /з) = 14.171190(рі -р4). В то же время отток жидкости из блока Ь2 равен -1.540083/?! + 5.661324/?2 - 2.581158/?3 - 1.540083/э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 Если М = AKÁr, то выполнены равенства 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 ■■

Г с

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 21228. 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 07632: 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 48855. 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 37999. 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.