Влияние численной диффузии на скорость роста вязких пальцев при численной реализации модели Писмана методом конечных объемов

Полный текст

Введение Процессы вытеснения более вязкой жидкости менее вязкой или газом широко распространены в современной нефтедобыче. Одна из технологий предполагает вытеснение нефти из пористого пласта смесью воды и полимера. При смешивающемся вытеснении менее вязкая вытесняющая жидкость прорывается через слой вытесняемой, образуя в ней неустойчивые структуры, называемые вязкими пальцами, или неустойчивостью Саффмана-Тейлора [20]. Несмотря на то, что точное решение задачи вытеснения с гладкой границей раздела существует, оно является неустойчивым, и любое возмущение на границе раздела сред приводит к росту «пальцев», что вызывает трудности точного описания параметров роста таких структур. Тем не менее в ряде случаев существует оценка длины вязких пальцев (см., например, [3]). Эта неустойчивость обычно исследуется не на примере фильтрационного процесса, а на примере течения в узком канале между двумя параллельными пластинами, в так называемой ячейке Хеле-Шоу. Это обусловлено как более простой организацией наблюдения при проведении экспериментов [5, 21], так и удобством построения возмущенных аналитических решений [14, 16-18], что упрощает верификацию и валидацию построенных численных моделей [10, 15, 17, 27]. Однако большую практическую значимость для нефтедобывающей отрасли имеет описание указанной неустойчивости именно в случае фильтрационного процесса [19, 22-24]. Несмотря на наличие применимой в ряде случаев аналитической оценки длины вязких пальцев (см., например, [20]), основным методом исследования остается численное моделирование [19, 22-24]. Как и в других задачах, одним из важнейших вопросов при использовании этого метода исследования является разделение двух типов эффектов: 1) присущих математической непрерывной модели явления и 2) обусловленных погрешностями аппроксимации уравнений и накоплением ошибки вычислений в ЭВМ. При численном описании процесса фильтрации сходным влиянием на скорость роста пальцев обладают модельная и численная диффузии. Модельная диффузия характеризуется совокупным действием жидкостной диффузии и дисперсии. При моделировании аппроксимация потока вносит в получаемое решение погрешности, влияние которых на решение похоже на воздействие физического диффузионного процесса. Такую погрешность можно определить как численную диффузию, которая накладывает определенные требования на выбор расчетной сетки. В работе представлены результаты численного моделирования вытеснения нефти смесью воды и полимера на основе модели Писмана [6, 7]. Численные эксперименты реализованы с помощью пакета DuMux, см. [13], представляющего собой программную библиотеку, предназначенную для моделирования нестационарных гидродинамических задач в пористых средах. Пакет программ использует вариант метода конечных объемов «vertex-centered». Исследовано влияние диффузии на скорость роста «вязких пальцев». Получены зависимости скорости переднего фронта (передней границы зоны смешивания) от значения модельной диффузии для трех моделей вязкости. Показано, что влияние численной диффузии на скорость роста «вязких пальцев» ставит ограничения на расчеты при малых значениях модельной диффузии. 16. Постановка задачи Рассмотрим участок пласта, заполненный нефтью. В начальный момент времени с левой границы рассматриваемого участка начинается закачка смеси воды и полимера как результат работы горизонтальной нагнетающей скважины, находящейся на большом расстоянии. Вытесняемая нефть уходит через правую границу. Поскольку уровень водонасыщенности в зоне полимерного заводнения мы считаем постоянным, соответствующим уровню поршневого вытеснения, то для описания процесса смешивающегося вытеснения используется традиционная однофазная модель Писмана: ⎧ ⎪⎨φ∂tc + div(u · c)= DΔc, u = -km(c)∇p, ⎪⎩div(u)= 0, (1) где φ - пористость породы, c - концентрация химического вещества в воде, u = u(u1, u2) - ско- 1 рость потока, D - коэффициент диффузии, k - проницаемость породы, m = μ(c) - функция обратной вязкости смеси, p - давление. Закон сохранения массы φ∂tc + div(u · c)= DΔc с учетом уравнения неразрывности div(u)=0 может быть записан следующим образом: φ∂tc + u∇c = DΔc. (2) Областью, в которой изучается решение системы уравнений (1), является бесконечная полоса Π(H)= {(x, y) ∈ R2 : y ∈ (0,H) ширины H. На горизонтальных границах полосы ставятся условия периодичности. Положим, что задано некоторое условие на постоянный поток при x = -∞: u(-∞, y)= q0e1, q0 ∈ R, где e1 - единичный орт по оси абсцисс. Тогда H r u1(t, x, y)dy = q0H ∀x ∈ R. (3) 0 Такое условие задает постоянную скорость закачки на нагнетающей скважине, которая в рассматриваемом случае предполагается расположенной на бесконечности. ВЛИЯНИЕ ЧИСЛЕННОЙ ДИФФУЗИИ НА СКОРОСТЬ РОСТА ВЯЗКИХ ПАЛЬЦЕВ В МОДЕЛИ ПИСМАНА 555 Заметим, что мы имеем следующую картину начала процесса смешивания: левая часть пря- 2 моугольника Π-(H) = {(x, y) ∈ R : x ∈ Π(H): x < 0 заполнена раствором воды и полимера с концентрацией cmin . Правая часть Π+ (H)= {(x, y) ∈ R2 : x ∈ Π(H): x > 0 заполнена нефтью с концентрацией cmax. В процессе заводнения средняя концентрация в направлении потока будет плавно меняться от cmin до cmax, при этом существенное отличие от крайних значений будет лишь в ограниченной, но растущей зоне смешивания (см. рис. 1). На схеме передний фронт Γf движется со скоростью vf , а задний фронт Γb - со скоростью vb. Рис. 1: Схема зоны смешивания. Черным цветом выделена зона, заполненная нефтью, серым - смесью воды и полимера. Зона смешивания ограничена пунктирными линиями Γf и Γb, которые соответственно являются передним и задним фронтами зоны смешивания. Fig. 1: Schema of the mixing zone. The zone lled with oil is highlighted in black, and the zone lled with a mixture of water and polymer is marked with gray. The mixing zone is limited by the dashed lines Γf and Γb,, which are respectively the leading and trailing fronts of the mixing zone. Для обезразмеривания использованы следующие характерные величины: x x˜ = , H y˜ = y , E = H q0 t, φH u u˜ = , q0 p˜ = k q0H p. (4) После обезразмеривания область Π (H) сведется к единичной полосе Π(1), скорость закачки q0 станет равна единице, а система уравнений (1)-(2) примет вид: 1 ⎧ ⎪∂˜c + ˜ (u˜ · c)= Δ˜ c, ⎨ t ∇ Pe u˜ = -m(c) ˜ p˜, (5) ∇ ⎪⎩ ˜ ∇u˜ = 0, где P e= q0H D - число Пекле [26]. 556 Д. Е. АПУШКИНСКАЯ, Г. Г. ЛАЗАРЕВА, В. А. ОКИШЕВ 17. Метод решения Численные эксперименты реализовались с помощью пакета DuMux, представляющего собой программную библиотеку, предназначенную для моделирования нестационарных гидродинамических задач в пористых средах. Пакет программ DuMux является открытым программным обеспечением и построен на оболочке DUNE (Distributed and Uni ed Numerics Environment), см. [13]. Данный пакет предназначен для эксплуатации только на операционных системах, входящих в семейство дистрибутивов Linux. В нашем случае использовалась Ubuntu 20.04, установленная на компьютер с четырехъядерным процессором Intel Core i5-10300H CPU @ 2,50 ГГц. Решение системы уравнений (5) основано на неявном методе конечных объемов «vertexcentered» и модели двухфазного потока 2pnc, см. [12]. Область расчета представляет собой прямоугольник длиной 60 м и шириной 31,415 м (высота пласта). Дискретизация осуществлялась с помощью прямоугольной сетки размером в 98766 ячеек, с шагом по времени τ =10 минут. Расчеты производились для трех моделей вязкости: § линейной: μ1(c)= A + B1c; § квадратичной: μ2(c)= A + B2c2; § экспоненциальной: μe(c)= A + Beekc; где A = 0,3 · 10-3, B1 = 3,8, B2 = 2533, Be = 0,017 · 10-3 и k = 3853. Вязкость рассчитывалась как функция концентрации полимера. Проницаемость пласта предполагалась изотропной. В целях стимуляции возникновения нестабильностей в динамической картине потока жидкости, приводящих к появлению вязких пальцев, была введена модельная вариация проницаемости пласта от 77,2 до 82,9 мД согласно закону логнормального распределения. Заданы следующие параметры для пласта, начальные и граничные условия: § пористость φ = 0,188; см3 § плотность нефти 0,773 г ; § вязкость нефти μn = 1,4 сП; • · на левой границе поставлено условие на постоянный поток раствора 3,68425 10-3 кг м · с , что сут соответствует площадной скорости закачки 10 м2 ; § начальное условие на водонасыщенность задается в соответствии с уровнем поршневого вытеснения и равна 0,5363. Для оптимизации процесса расчета по времени параллельно запускалось четыре симуляции, каждая в своем терминале, т. е. на одну симуляцию выделялось одно ядро процессора. Это было связано с ограниченными ресурсами компьютера и необходимостью провести наибольшее количество численных экспериментов. Надо отметить, что помимо высоких требований к вычислительной мощности компьютера, потребовался большой объем памяти для хранения полученных данных. При запуске симуляции с помощью MPIrun, время расчета сократилось не более чем в полтора раза. Таким образом, параллельный запуск на четырех ядрах работал быстрее, чем последовательный запуск четырех симуляций с использованием MPIrun. Однако запуск MPIrun не принес существенного преимущества в скорости вычислений. 18. Результаты Для проведения расчетов в конфигурационных файлах WettingPhaseViscosityModel и PolymerDi usionCe cient изменялись модель вязкости и значение модельной диффузии. Визуализация реализована с использованием приложения Paraview. С ростом модельной диффузии вязкие пальцы размываются (см. рис. 2, 3) и при значениях D 10-4 пропадают. В этом случае мы имеем поршневое вытеснение (см. рис. 2 d). На рис. 4 показано, что длина пальцев в один и тот же момент времени при D = 10-9 значительно больше, чем при D = 10-6. Также можно отметить, что скорость переднего фронта постепенно увеличивается с ростом модельной диффузии. Для определения скорости роста пальцев необходимо определить зависимость положения переднего и заднего фронтов от времени (см. рис. 5 a, b, c, d). Скорость переднего фронта можно оценить как коэффициент наклона линейного приближения соответствующих графиков (cм. рис. 5). Заметим, что при значениях модельной диффузии D 10-4 график зависимости ВЛИЯНИЕ ЧИСЛЕННОЙ ДИФФУЗИИ НА СКОРОСТЬ РОСТА ВЯЗКИХ ПАЛЬЦЕВ В МОДЕЛИ ПИСМАНА 557 (a) D =0 (b) D = 10-9 (c) D = 10-6 (d) D = 10-4 Рис. 2: Карты концентраций для различных значений модельной диффузии D при t = 626525 c (экспоненциальная модель вязкости). Fig. 2: Concentration maps for di erent values of model di usion D for t = 626525 sec (exponential viscosity model). (a) D =0 (b) D = 10-6 (c) D =0 (d) D = 10-6 Fig. 3: Карты концентраций для линейной (a, b) и квадратичной (c, d) моделей вязкости при различных значениях модельной диффузии D при t = 907350 c. Fig. 3: Concentration maps for linear (a, b) and quadratic (c, d) viscosity models for di erent values of model di usion D for t = 907350 sec. 558 Д. Е. АПУШКИНСКАЯ, Г. Г. ЛАЗАРЕВА, В. А. ОКИШЕВ Рис. 4: Наложение карт концентраций при D = 10-9 и D = 10-6 при t = 626525 c (экспоненциальная модель вязкости). Fig. 4: Overlay of concentration maps for D = 10-9 and D = 10-6 for t = 626525 sec (exponential viscosity model). ⎫ ⎫ ⎬ ⎬ ⎭ ⎭ ⎫ ⎫ ⎬ ⎬ ⎭ ⎭ (a) (b) ⎫ ⎬ ⎭ ⎫ ⎬ ⎭ (c) (d) (e) (f) Fig. 5: Зависимость положения фронтов от времени: для линейной (a), квадратичной (b) и экспоненциальной (c) функций вязкости, для D = 10-4(d), переднего фронта для D = 0 (e) и для D = 10-6 (f). Fig. 5: Dependence of front positions on time: for linear (a), quadratic (b), and exponential (c) viscosity functions, for D = 10-4(d), leading front for D =0 (e) and D = 10-6 (f). ВЛИЯНИЕ ЧИСЛЕННОЙ ДИФФУЗИИ НА СКОРОСТЬ РОСТА ВЯЗКИХ ПАЛЬЦЕВ В МОДЕЛИ ПИСМАНА 559 положения заднего фронта от времени становится параллельным оси абсцисс (см. рис. 5 d), что говорит нам об увеличении длины вязких пальцев. С увеличением модельной диффузии скорости фронтов постепенно возрастают для любой модели вязкости. В случае нулевой модельной диффузии (см. рис. 5 e) положения передних фронтов сильно отличаются для разныхмоделей вязкости. Причем для экспоненциальной модели вязкости скорость переднего фронта выше, а для линейной ниже, чем для квадратичной. Углы наклона графиков, а значит, и скорости движения фронтов, сильно различаются. Аналогичным образом соотносятся положения фронтов для модельной диффузии D = 10-6 для разных моделей вязкости (см. рис. 5 f). При этом графики положения передних фронтов для линейной и экспоненциальной моделей вязкостей параллельны, а для квадратичной и экспоненциальной совпадают в начале расчета. Это показывает, что скорости передних фронтов при D = 10-6 почти совпадают для всех трех моделей вязкости. Графики зависимости скорости переднего фронта от значения модельной диффузии для разных моделей вязкости приведены на рис. 6. Они показывают, что рост скорости переднего фронта при достаточно малом коэффициенте диффузии нарушает монотонный характер зависимости скорости от диффузии. Количественные различия скоростей составляют величины не более единицы (см. рис. 6 b). Более того, минимум скоростей достигается при D = 10-6-10-5 во всех случаях. Возможно, этот факт говорит о том, что численная диффузия больше этой величины. Поэтому добавка модельной диффузии не вносит изменений в результат вычислений. Вычисления на более грубых сетках показали, что эта тенденция сохраняется. 1. (b) Рис. 6: Зависимость скорости переднего фронта от диффузии в диапазоне от 0 до 0,1 (a) и от 0 до 10-4 (b). Fig. 6: Dependence of the leading front velocity on di usion in the range from 0 to 0.1 (a) and from 0 to 10-4 (b). Грубо оценить численную диффузию Dnum можно по формуле: Dnum = u Δx . Такая априор- 2 ная оценка согласуется с наблюдаемыми результатами, но требует уточнений для используемого метода конечных объемов «vertex-centered». В ряде случаев существуют более тонкие оценки величины численной диффузии. Например, известна оценка влияния численной диффузии на решение уравнений Навье-Стокса при моделировании поверхностных волн в случае применения метода конечных объемов. В работе [4] предлагается методика оценки численной диффузии, выражаемой коэффициентом уменьшения амплитуды волны при прохождении ею одной своей длины (коэффициентом затухания). Вывод более тонкой оценки численной диффузии для решения системы (5) является следующим необходимым шагом в дальнейших исследованиях. 19. Выводы Рассмотрена численная модель вытеснения нефти смесью воды и полимера на основе модели Писмана. Проведены численные эксперименты с помощью пакета DuMux, представляющего собой программную библиотеку, предназначенную для моделирования нестационарных гидродинамических задач в пористых средах. Пакет программ использует вариант метода конечных 560 Д. Е. АПУШКИНСКАЯ, Г. Г. ЛАЗАРЕВА, В. А. ОКИШЕВ объемов «vertex-centered». Исследовано влияние диффузии на скорость роста «вязких пальцев». Получены зависимости скорости переднего фронта от значения модельной диффузии для трех моделей вязкости. Показано, что влияние численной диффузии на скорость роста «вязких пальцев» ставит ограничения на расчеты при малых значениях модельной диффузии.

Об авторах

Д. Е. Апушкинская

Российский университет дружбы народов

Автор, ответственный за переписку.
Email: apuskinskaya-de@rudn.ru
Москва, Россия

Г. Г. Лазарева

Российский университет дружбы народов

Email: lazareva-gg@rudn.ru
Москва, Россия

В. А. Окишев

Российский университет дружбы народов

Email: okishev-va@rudn.ru
Москва, Россия

Список литературы