Application of the finite element grouping procedure to improve the efficiency of unsteady multiphase flow simulation in.pdf Современные программные комплексы численного ЗО-моделирования процессов многофазной фильтрации в пористой среде широко используются при разработке нефтегазовых месторождений. Сложное строение коллекторов требует создания эффективных вычислительных схем, позволяющих проводить расчеты для моделей реальных месторождений с большим числом слоев, сильной неоднородностью среды и с большим количеством работающих скважин и зон перфораций. Успешность, а иногда и сама возможность решения таких важных для практики задач, как построение гидродинамических моделей месторождений путем решения обратных задач [1, 2] и синтез оптимального управления разработкой [3, 4], во многом зависит от эффективности вычислительных схем, применяемых при решении прямых задач. Часто для решения задач фильтрации используют конечно-разностные численные схемы. В частности, подобные подходы используются в широко распространенных коммерческих симуляторах нефтедобычи, таких как Eclipse, Tempest и др. Однако многие авторы отмечают, что методы конечных разностей и конечных объемов обладают недостаточной геометрической гибкостью и недостаточной точностью при моделировании многофазных течений в высоконеоднородных средах [5, 6]. Публикуется довольно много работ, в которых для решения задач в сложных высоконеоднородных средах предлагается использовать различные модификации метода конечных элементов (МКЭ), например совместное использованием МКЭ и метода конечных объемов (FEFVM) [7, 8], mixed FEM и его специальные модификации (mixed hybrid FEM) [9-11]. Однако при использовании таких методов вычислительные затраты при решении ориентированных на практику задач нефтедобычи со сложной (многослойной латерально неоднородной) структурой среды и большим числом действующих скважин могут быть очень велики. Проблема снижения вычислительных затрат особенно актуальна при решении обратных задач, когда в процессе восстановления модели сложного месторождения требуется решать большое число соответствующих прямых задач. В данной работе будет использован вариант МКЭ с непрерывными базисными функциями Гильбертова пространства функций с производными, суммируемыми с квадратом (Continuous Ga-34 Применение процедуры группирования конечных элементов lerkin [2, 12]). Необходимая точность выполнения законов сохранения масс достигается за счет специальной процедуры балансировки потоков [13]. Данная вычислительная схема была верифицирована на тесте SPE-10 [14] и применялась при моделировании реальных месторождений высоковязкой нефти Республики Татарстан [2]. Повышение вычислительной эффективности для задач такого класса достигается в том числе и за счет использования специальных некомфорных сеток [15], которые позволяют существенно сократить число степеней свободы в конечноэлементных аппроксимациях без увеличения погрешности численного решения. Свою эффективность предложенный подход показал при решении обратных задач [2] и при синтезе оптимального управления разработкой месторождения [4]. Принцип построения предлагаемой вычислительной схемы аналогичен принципу известных схем IMPES [16] и IMPEC [10, 11] с конечноэлементным (неявным) расчетом давления и явным переносом фаз между ячейками конечноэлементной сетки на каждом временном шаге. Вычислительные затраты напрямую зависят от величины временного шага, дробление которого может потребоваться при большой неоднородности среды и сильно изменяющейся в пространстве скорости потока. В данной работе будут представлены схема группирования конечных элементов по временному шагу и алгоритм моделирования перетоков фаз с учетом разбиения ячеек сетки на группы, позволяющие минимизировать количество вычислений поля давления, не понижая качества аппроксимации многофазного потока. В качестве теста использовалась задача сравнительного проекта SPE-10. Показана сходимость данного метода, проведены сравнение с данными участников проекта SPE-10, анализ затрат машинного времени и точности получаемых решений при разных временных шагах. 1. Математическая модель В задачах моделирования месторождений расчет движения многофазной смеси выполняется в неоднородной пористой среде, которая характеризуется зависящими от пространственных координат тензором структурной проницаемости К и пористостью Ф. Скорость ѵт движения в пористой среде каждой фазы фильтрующейся смеси подчиняется закону Дарси: т / j \\ Г” =-b-K(^grad(p+p;”)+(o,o,p"g) J, где т - номер фазы, ѵ1" - скорость потока фазы, Р - давление, Р'" - капиллярное давление, g -ускорение свободного падения, к™ - коэффициент относительной фазовой проницаемости, цт - динамическая вязкость, рт - плотность фазы. Распределение давления в расчетной области Q описывается краевой задачей [2]: ( -div I кт К (grad (P + PC ) + ( 0,0, рт g )' и=1 л V ' ѵ ' P Г1 =рг, т=1 Л' Iк™ К (grad (P + Рт) + (0,0, pmg)‘ т=1Л V Ѵ Ѵ NP _ = ІГ’°, т=1 NP ■n=ZfmX- Г2 т=1 (1) (2) (3) Здесь NP - количество фаз. Функции f соответствуют объемным источникам (стокам) фаз в расчетной области Q, возникающим, например, вследствие выделения / поглощения газа жидкими фазами, химических процессов или сжатия / расширения фаз при изменении давления. Г1 - объединение тех границ расчетной области Q, где задано давление Рг, а Г2 - объединение границ Q, где Г NP задан поток смеси f = I /т’ . Функция f1 не равна нулю на тех границах из Г2, которые соотт=1 ветствуют активным (в интервале времени At) зонам перфорации. Остальные границы из Г2 являются непроницаемыми (на них fГ = 0). 35 М.Г. Персова, Ю.Г. Соловейчик, И.И. Патрушев, А.С. Овчинникова Краевая задача (1)-(3) решается методом конечных элементов на шестигранных неконформных сетках [15]. По полученным значениям давления вычисляются объемы смеси, перетекающие через грани Г; конечных элементов Qe за единицу времени: Qt = -1 Г,. 2^ к( grad (P + Pcm ) + (о,0, pmg )T m=1 p V где йг q - внешняя (по отношению к Q„) нормаль к Г;. Для внутренних граней Г; объем смеси, перетекающий за единицу времени из конечного элемента (ячейки) в смежную с ним по грани Г; ячейку Q.k (Г; =fien£\\) определяется как взвешенное среднее объемов, перетекающих через эту границу на конечных элементах Q.e и Qk : Qt i Ь к Ьк + Ье Ь е Ьк +Ье Qr (4) ~ Л(Р yjn ^ Коэффициенты А, в (4) ().к для Q,. и /.„ для С1е) определяются как Х = КУ]--, где К = Щ. КЯГ , т=1 Г| " " т.е. К определяется по значению тензора К на соответствующем конечном элементе. Если для решения краевой задачи (1)-(3) применяется МКЭ с базисными функциями из H1 (так называемый CG [12]), то получается численное решение, не гарантирующее сохранения масс веществ в фильтрующейся смеси [7, 17] (закон сохранения в этом случае лишь аппроксимируется с той или иной точностью в зависимости от подробности сетки). Поэтому мы используем специальный метод балансировки потоков [13], который корректирует перетекающие объемы Qt таким образом, чтобы законы сохранения масс отдельных фаз (и их компонент) были выполнены с необходимой точностью. По сбалансированным потокам смеси QT осуществляются перетоки фаз через Г; и вычисляется новый фазовый состав в конечных элементах. Для этого определяются объемы фаз , перетекающих через грань Г в единицу времени. В условиях, когда эффект гравитации или капиллярного давления является существенным, разные фазы могут перетекать через грань Г в противоположных направлениях. С учетом этого алгоритм выполнения перетоков между конечными элементами выглядит следующим образом. Для всех конечных элементов и принадлежащих им граней Г вычисляются численные по токи фаз Qm, Q по формуле Q m Г Q m -1 - K j^grad (P + Pcm) + (0 0 Pmg) ' ”r,,Qe . (5) Уравнение (5) определяет не только величину потока m-й фазы через грань Г , но и его направление по отношению к конечному элементу Qe. Фаза m вытекает из конечного элемента Qe, если значение Qm q положительное, и втекает при отрицательном значении. Численный поток фазы Qm через грань Гі берется с подветренной стороны: Qmt = Г,Qe если QTi, Qe > 0, (6) Qm,^ если Qm,^ < 0. Затем для граней Г; для каждой фазы вычисляются величины Dm,, которые фактически определяют долю m-й фазы в потоке смеси через грань Г : э (7) om = Qmi/itei / n=1 36 Применение процедуры группирования конечных элементов Далее мы корректируем численные потоки фаз Q™ таким образом, чтобы сумма потоков всех фаз была строго равна сбалансированному потоку смеси QT , перетекающему через грань Г;. Для этого вычислим разницу между численным и сбалансированным потоками смеси и распределим ее между потоками фаз пропорционально величинам D™ , т.е. перевычислим Qm по формуле ( ~ NP \\ 0£=0£+\\0г-Ш • Dm (8) Таким образом, объем m-й фазы V™, который за время At перетекает через грань Г;, вычисляется с использованием по формуле Vrm = Qmi At, (9) и новые значения насыщенностей на каждом элементе Qe на конец шага по времени At о т _ - Sq mes(Oe) + AF£ + Z V™ - Z Цг" /(mes(Oe)), (10) ielm’m ieiQ, где mes (Qe) - объем ячейки Qe; AVm - дефицит / профицит объема m-й фазы в конечном элементе Qe, который может образоваться, например, в результате химических процессов или сжатия / расширения фаз при изменении давления; Iq1’™ , Iq™ - множества номеров граней элемента Qe, через которые m-я фаза вытекает из Q и, соответственно, втекает в Q. На основе полученных значений вычисляются новые значения фазовых проницаемостей кт согласно заданным зависимостям к™ от насыщенностей фаз. Затем осуществляется переход к следующему шагу по времени, на котором процедура повторяется, начиная с расчета давления. Теперь более подробно остановимся на выборе временного шага At. Кроме того, что он непосредственно влияет на точность аппроксимации по времени, его значение должно быть ограничено величиной объема фаз в ячейках, из которых эти фазы вытекают. Шаг At должен быть таким, чтобы в каждой ячейке суммарный вытекающий объем фазы m не превышал имеющийся объем подвижной фазы в ней. Это естественное условие определяет его предельное (максимальное) значение: At < (( г 1 в 2s-1 раза меньше Atmain . Количество таких групп обозначим NG . Для того, чтобы распределить ячейки Qe по группам, для каждой m-й фазы в ячейке Qe согласно условию (11) определяем временной шаг, являющийся допустимым для этой фазы: ) Z Qmt iel tom* = (К -S"Ce)mes(П.)Ф + Д^г:) Z Qm. (12) Номер группы для конечного элемента Qe выбирается из условия (13) minjg : ЛУ ^ min Atm ,n. I . s m ' Перенос фаз между ячейками осуществляется следующим образом. Обработка ячеек начинается с группы с самым большим номером g = NG (т.е. с самым маленьким шагом по времени). Делается два шага по времени Atr для всех ячеек этой группы. Затем делается один шаг по времени для ячеек группы с номером g = NG -1, после чего опять делается два шага по времени для ячеек группы с номером Ng . Далее делается еще один шаг по времени для ячеек с номером группы g = NG -1. После этой второй обработки данной группы делается первый шаг обработки группы с номером g = NG - 2 и опять осуществляется возврат к обработке группы с номером g = NG . Начало Цикл по ячейкам Qe из группы G (е е I ) ►I Цикл по граням Г (г е Іп ) ► CurGf^l - CurG|g] +1 г = номер группы соседнею элемента Qk (Qk nQe = 12 ) Цикл по граням Г ячеек Qe из группы Gg (е е / , і е ) Нет г = номер группы соседнего элемента П, (Ц,. n Q(, = Гі ) I Іет Расчет насыщенностей Конец Цикла по гранямг/] Вычисление потоков Конец Конец Цикла по граням Г ячеек L2 из группы & Выбор нового номера группы (ур. 12-13) Вычисление ооъемов Нет CurGfgl = О Конец Цикла по ячейкам і\\ из группы и . Рис. 1. Алгоритм, реализующий расчет перетоков фаз Fig. 1. Algorithm of phase transport 38 Применение процедуры группирования конечных элементов Таким образом, после первой обработки каждой группы с номером g осуществляется возврат к обработке ячеек группы с номером NG (и далее до группы g), а после второй обработки группы с номером g осуществляется переход к группе g - 1 (с обнулением счетчиков CurG обработки групп с большими номерами). Алгоритм завершается после обработки группы с номером 1 (т.е. группы с самым большим шагом AtG ).Обработка ситуации, когда фазы перетекают между ячейками из разных групп, осуществляется следующим образом. Прежде всего отметим, что при вычислении перетекающих объемов фаз Ур по формуле (9) всегда используется шаг Atr группы того элемента, из которого вытекает поток QP . При пересчете насыщенностей в ячейках это позволяет использовать вытекающие из них объемы фаз Ур без изменений. А втекающие объемы вычисляются в зависимости от того, какой группе принадлежит соседний элемент, из которого фазы перетекают в обрабатываемую ячейку. Таким образом, если при обработке некоторого конечного элемента из группы Gg через одну из его граней втекала фаза из ячейки группы с меньшим или равным номером r , то перетекающий объем УР делится на 2g-r, чтобы учесть разницу во временных шагах. Если же в элемент из группы Gg втекает фаза из ячейки группы с большим номером r , то используется ранее накопленное значение перетекающего объема УттР , которое было получено в ходе предыдущих обработок 1 i ячейки группы с номером r . Алгоритм, реализующий процедуру перетоков смеси для расчета нового состояния ячеек с учетом группирования, изображен на рис. 1. 3. Численные эксперименты Продемонстрируем работоспособность и эффективность процедуры группирования на примере десятой тестовой задачи из сравнительного проекта SPE, рассмотренной в [14]. Модель представляет собой куб с размерами 1200 х 2200 х 170 фут3 (365,76 х 670,56 х 51,816 м3), в углах которого расположены четыре добывающие скважины (P1-P4), а в центре находится нагнетательная скважина. Добывающие скважины работают при заданном в них давлении 272 атмосферы, а нагнетательная скважина закачивает 5 000 баррелей воды в сутки. Исходная модель содержит 85 слоев и характеризуется высокой неоднородностью структурной проницаемости и пористости. Значения параметров задачи заданы на регулярной сетке, содержащей 1 122 000 ячеек (60 х 220 х 85 ячеек), и варьируют в диапазоне от 6,65 • 10-4 до 20 000 мД для коэффициентов тензора проницаемости и в диапазоне от 1,4 10-6 до 0 для пористости. Для этой модели была выполнена предварительная обработка с укрупнением ячеек сетки по вертикали (upscaling). В результате расчет проводился на сетке из 244 052 ячеек. Время жизни «месторождения» - 2 000 суток. Для оценки эффективности предлагаемого нами подхода с группированием ячеек была проведена серия расчетов с разными временными шагами Atmam, равными 100, 50, 20, 10, 5, и 1 сут. Также был проведен расчет без использования процедуры группирования. В нем расчет давления выполнялся с шагом At = 20 суток, а перетоки фаз - с одинаковым для всех ячеек временным шагом, удовлетворяющим условию (11). На рис. 2 представлены соответствующие результаты для скважин P1 и P3 (рис. 3): дебит нефти в баррелях (см. рис. 2, а) и обводненность (см. рис. 2, b) для шагов Atmam, равных 100, 50, 20 и 1 сут. Очевидно, с дроблением шага по времени наблюдается сходимость. Данные дебита нефти на скважине P3 при шаге Atmam на шагах 20, 10 и 5 сут отличаются от значений при A tmam=1 сутки в среднем на 2, 1 и 0,5% соответственно, а на скважине P1 - 1, 0,5 и 0,2%. Также от-39 М.Г. Персова, Ю.Г. Соловейчик, И.И. Патрушев, А.С. Овчинникова метим, что дебит нефти при расчете с шагом Afmam=1 сут очень хорошо совпадал с результатами, полученными без применения процедуры группирования. Аналогичная картина сходимости наблюдается и для обводненности добываемой смеси. Рис. 2. Графики дебита нефти (а) и обводненности (b) для расчетов с различными At” : квадрат - 100 сут, незакрашенный круг - 50 суток, треугольник - 20 суток, закрашенный круг - 1 сутки Fig. 2. Oil rate (a) and water cut (b) curves for calculation with different At™” : square - 100 days, open circle - 50 days, triangle - 20 days, filled circle - 1 day a) b) Рис. 3. Распределение ячеек по группам (зеленый цвет - группа с большим номером, желтый цвет - с меньшим) Fig. 3. Distribution of cells into groups (green is the group with a higher number, yellow is the group with a lower number) Временные затраты для расчетов с различным шагом ytman Atma'n, сут Количество групп Расчет давления Расчет потоков Расчет перетоков Общее время расчета 1 8 13 ч 56 мин 1 ч 55 мин 1 ч 14 мин 17 ч 37 мин 5 10 2 ч 40 мин 24 мин 31 мин 3 ч 42 мин 10 11 1 ч 25 мин 13 мин 29 мин 2 ч 12 мин 20 12 43 мин 7 мин 27 мин 1 ч 20 мин 50 14 17 мин 3 мин 26 мин 48 мин 100 15 8 мин 2 мин 26 мин 38 мин Без группирования 43 мин 6 мин 66 ч 46 мин 84 ч 31 мин 40 Применение процедуры группирования конечных элементов В таблице приведены вычислительные затраты основных процедур описанной вычислительной схемы. На каждом временном шаге Atmain производится решение двух СЛАУ: первой - при решении краевой задачи (1)-(3) методом конечных элементов (размер СЛАУ определяется количеством узлов сетки), второй - в ходе выполнения процедуры балансировки для расчета потоков смеси (размер этой СЛАУ определяется количеством граней [13]). Для решения обоих СЛАУ применялся прямой решатель PARDISO из библиотеки Intel MKL. Сопоставляя погрешность по времени с вычислительными затратами при различных значениях Atmam, можно считать оптимальным значение Atmain = 20 суток. Для этого расчета на рис. 3 приведено распределение ячеек по группам на конец расчета, где желтый цвет соответствует группе с наименьшим номером (шаг по времени AtG = 20 сут), а темно-зеленый цвет соответствует ячейкам из группы с наибольшим номером (шаг по времени Atr = 28 мин). Из рисунка видно, что ячейки, расположенные вдоль основных направлений течения (от нагнетательной скважины I1 к добывающим P1, P2, P3, P4) и в зонах повышенной проницаемости, попадают в группы с большими номерами, и основные вычислительные затраты при расчете перетоков приходятся на них. Если же не использовать процедуру группирования, то затраты машинного времени составляют порядка 85 час из-за условия Куранта-Фридрихса-Леви (CFL) и связанной с ним необходимости задания довольно мелкого шага по времени для выполнения процедуры перетоков фаз (At примерно 20-30 мин). Заключение Рассмотрена вычислительная схема моделирования многофазных потоков в пористых средах при решении задач нефтедобычи. Схема основана на неявном расчете давления с использованием МКЭ, балансировке численных потоков смеси и явном переносе фаз между конечными элементами. Предложена процедура группирования и основанный на ней алгоритм переноса фаз, позволяющие выполнять пересчет насыщенностей в ячейках с разными шагами по времени в зависимости от скорости течения и объема фаз в них. Тем самым устраняется необходимость использования слишком мелких шагов по времени из-за возможного нарушения условия CFL на отдельных ячейках, и резко снижаются вычислительные затраты без ущерба точности получаемого решения. Эффективность предложенной процедуры группирования и алгоритма переноса фаз продемонстрирована на тесте SPE-10, в котором нефтедобыча моделируется в высоконеоднородной среде. Для этой задачи проведенные вычислительные эксперименты показали, что за счет использования предложенной процедуры группирования время счета сокращается примерно на полтора порядка без снижения точности получаемого решения.

Скачать электронную версию публикации