УДК 519.633 + 004.021 +532.546
А. И. Шангараева
АНАЛИЗ АЛГОРИТМА УСКОРЕНИЯ РАСЧЕТА НЕФТЕНАСЫЩЕННОСТИ
В ОДНОМЕРНОМ СЛУЧАЕ
Ключевые слова: вычислительный метод, ускорение гидродинамических расчетов, компьютерное моделирование, число
Куранта.
Данная статья освещает вопрос создания и обоснования алгоритма для ускорения численных расчетов параболических уравнений на примере уравнения для нефтенасыщенности.
Keywords: numerical algorithms, the acceleration of hydrodynamic computations, computer simulation, Courant number.
This article is about the algorithm, which allows to increase the speed of calculations ofparabolic equations.
В данной статье предлагается метод ускорения численных расчетов уравнений переноса по явной схеме с использованием локального числа Куранта. Для подтверждения достоверности результатов приводится сравнение аналитического и численного (с ускорением и без ускорения) решений в одномерном случае.
Рассмотрим случай изотермической фильтрации двухфазной жидкости в пласте, вскрытом системой добывающих и нагнетательных скважин. Для простоты изложения алгоритма рассмотрим случай отсутствия капиллярных сил и сил тяжести. Фильтрация жидкостей происходит по обобщенному линейному закону Дарси. Модель описывается системой [1]:
pJi _v(k .k'-vp ) = о
ds
dp
(1)
m— = -V(f -u)-~Pw dt ' w dt
где: 5 - насыщенность, р - давление, ¡3 -
суммарная упругоемкость среды, 3к - эффективная упругоемкость системы, относящаяся к водной фазе, к - абсолютная проницаемость, к ■ к - общая проводимость двухфазной системы, т -пористость, и = -к ■ к * ■ Ур - общая скорость флюидов, / = / (5, х) - функция Баклея-Леверетта
или доля воды в потоке.
Известно [2, 3], что уравнения данного типа на практике широко используются при гидродинамическом моделировании нефтяных месторождений. При этом множество различных параметров, влияющих на процедуру адаптации, создает необходимость множества расчетов. Временные рамки для выбора подходящей модели требуют ускорения компьютерных расчетов.
В случае 1МРБ8 [2, 3] технологий первый этап проводится по неявной схеме (расчет поля давления и скоростей), а для расчета второго этапа используется явная схема (расчет поля нефтенасыщенности). Такая технология
относительно проста в реализации и легко адаптируется под любые физические модели. Однако при использовании явной схемы возникает ограничение на шаг по времени в виде условия Куранта [1] (С < 1). Кроме того, ряд необходимых
вычислительных процедур (измельчение сеток и т.п.) приводит к еще меньшим временным шагам. Во многих случаях область около скважин рассчитывается специальным образом [2, 4].
Специфика задач моделирования
гидромеханики нефтяных месторождений состоит в том, что есть области с очень высокими и очень малыми скоростями, это, соответственно, потоки вблизи скважин и вдали от них. Условие Куранта может быть очень жестким за счет потоков в окрестности скважин, что повлечет за собой большие затраты по времени для расчета.
Идея состоит в усовершенствовании алгоритма явной схемы. В работе исследуется использование разных шагов по времени в разных областях месторождения так, чтобы в каждой области шаг был: во-первых, достаточно велик для уменьшения расчетных процедур и, во-вторых, удовлетворял условию Куранта для соответствующих ячеек. Подобная процедура впервые была изложена в работах П. П. Матуса [5], где предложено использовать кратный измельченный шаг по времени в области с высокими градиентами.
В [6] предложено обобщить подход на случай нескольких взаимосвязанных областей, требующих своего коэффициента измельчения. Целью данной работы является верификация алгоритма в одномерном случае и анализ его эффективности в радиальном случае.
Классический подход, в случае использования явной схемы, предполагает выбор шага по времени, равного минимальному среди всех шагов для каждой ячейки. Для начала вводится число Куранта [1] на основе численного аналога второго уравнения системы (1):
Vimi
gn+1 — gn
= —f iu iS i +
+f 1u 1S 1 — Pw
i +— i+— i +— 2 2 2
— pN
-V
которое равно:
u г
C = F——
m■ Vi
где V , £ - объем и боковая поверхность расчетной ячейки, Г - множитель запаса, равный
т
i— 1— 1— 2 2 2
максимальному значению /' ), индексом /
обозначены величины, относящиеся к конкретной I -й ячейке.
Наилучший вариант - для каждой ячейки иметь свой шаг по времени, удовлетворяющий условию Куранта С < 1:
СтУ.
Fu,
Однако непосредственное задание шагов по времени в каждой ячейке невозможно из-за необходимости согласованного расчета поля насыщенности. Вследствие этого предлагается метод, позволяющий использовать различные временные шаги для разных подобластей с их динамическим выбором по полю скорости фильтрации. Для каждой ячейки выбирается шаг по времени, равный минимальному из шагов, удовлетворяющих условию Куранта, соответствующих потокам, протекающим через грани ячеек (рис. 1):
{Vm )
F
-1/2;и 1+1/2
X
Рис. 1 - Иллюстрация процесса пересчета
1. Выбирается наименьший по всем ячейкам шаг по времени
гтт = ттг,.
1
2. Для всех ячеек формируется коэффициент увеличения шага
г.
тт
(2)
3. Коэффициенты для ячеек округляются в сторону нуля до числа, составляющего целую степень числа два (1, 2, 4, ...), тем самым определяется количество мелких шагов, на которых будет проведен один перерасчет для ячейки:
NSt = 2[log2 Ki 1 (3)
4. В каждой ячейке вычисляется свой шаг по времени
it= NSтт1П
5. Значение функции Баклея-Леверетта на каждой границе f.+1/2 пересчитывается один
раз за max (INS.; NSj+1) шагов.
Пункт 5 алгоритма объясняется следующим. Во-первых, при этом потоки переносимой величины будут согласованы, выполнено условие консервативности схемы. Во-вторых, условие устойчивости будет выполнено.
Поясним оба утверждения на примере. Пусть, например, Ti_1 = 2тт1п и т = 8тт1П, т.е. за 8 шагов т :
min
а) i _ 1 -ая ячейка будет пересчитываться 4 раза (1 раз в 2 шага) с шагом 2тт1П ;
б) -ая ячейка будет пересчитываться 1 раз с шагом 8тт1П . Соответственно, из (i _ 1) -ой ячейки
«вытечет» вправо (fu) • 2тт1П • 4, а в ячейку i
«втечет слева»
{fu ),-_L ' 8rmm и только потом будет
пересчитан.
Осталось показать, что пересчет раз в 8 шагов не приводит к потере устойчивости. Это следует из того, что шаг 8гтш в / -ой ячейке был выбран как минимальный, в том числе с учетом потока из (/ -1) -ой ячейки.
Проведено также исследование возможности выбора разных значений Е в разных ячейках. В результате выявлено: для выбора локального Е , заранее должна производиться оценка продвижения фронта, что является достаточно сложной задачей. Численные эксперименты показали, что определение для каждой ячейки не имеет смысла. Актуальным является лишь учет разных потоков и объемов ячеек.
Проведено сравнение вышеописанного варианта численного расчета с аналитическим решением для плоского случая [7]. Результаты расчетов достаточно хорошо совпадают с аналитическими результатами (рис.2).
Рис. 2 - Результат реализации аналитического и численного решений
Алгоритм проверен так же на примере задачи движения фронта насыщенности в радиально-симметричном случае (рис.3). Данные графики иллюстрируют совпадение расчетов без ускорения и с ускорением. Очевидно также, что движение фронта в обоих случаях идентично.
Покажем, как в данном примере выглядит коэффициент увеличения шага (рис.4).
Сплошной линией показан математический идеал, пунктирной линией - идеал реализации, полученные по формулам (2) и (3), соответственно.
2
Рис. 3 - Движение фронта насыщенности через равные шаги по времени в радиальном случае (сплошная линия - результат без ускорения, пунктирная линия - с ускорением)
£ . i "
Рис. 4 - Коэффициент увеличения шага по времени в зависимости от шага (сплошная линия - математический идеал, пунктирная линия -идеал реализации)
300
VV1
□
qj
7Г>П
4J
Ш
О ISO
H
tl ion
Ц
tn 50
0
Рис. 5 - Количество расчетов в зависимости от радиуса. Пунктирная линия - без ускорения, сплошная линия - количество пересчетов с ускорением. Точками показана кривая математического идеала
На рисунке 5 показано как уменьшается количество пересчетов при использовании
алгоритма ускорения. Из графика видно, что суммарное количество расчетов меньше примерно в 4 раза, то есть можно ожидать четырехкратного ускорения.
Численные расчеты показали выигрыш реального (физического) времени расчетов нефтенасыщенности примерно в 1,7 раза. Отличие реального времени от прогнозируемого объясняется тем, что в предложенной реализации специальных процедур оптимизации алгоритма не проводилось. Естественно предположить, что такая оптимизация позволит приблизиться к ожидаемому ускорению. Кроме того, по предварительным оценкам, трехмерный случай должен дать дополнительный выигрыш по времени.
Итак, предложен алгоритм, позволяющий уменьшить время расчета нефтенасыщенности в нефтяном пласте. Корректность работы алгоритма доказана проверкой в одномерном плоскопараллельном и радиально-симметричном случаях. Обоснована возможность ускорения в 1,7 раза, а потенциальная возможность - в 4 раза.
Литература
1. Х. Азиз, Э. Сеттари Математическое моделирование пластовых систем. М., Недра, 1982. - 407 с.
2. VIP-EXECUTIVE technical reference guide - режим доступа http://esd.halliburton.com/support/ LSM/ResMgmt/NexusVIPDT/Nexus/5000/5000_4/Help/Te chref.pdf.
3. Р. Д. Каневская Математическое моделирование гидродинамических процессов разработки месторождений углеводородов - Москва-Ижевск: Институт компьютерных исследований, 2002, - 140 с.
4. Д.В. Шевченко, О.В. Панченко Учет в численных алгоритмах последствий применения методов значительного улучшения свойств призабойной зоны скважин. // Вестник казанского технологического университета. 2014г., Т.17, №6. - С. 267-269.
5. П.П. Матус Консервативные разностные схемы для параболических уравнений второго порядка в подобластях. // Дифференциальные уравнения. 1993г., Т. 29, № 4. - С. 700-710.
6. Шевченко Д. В., Шангараева А. И., Использование локального числа Куранта для ускорения расчета нефтенасыщенности по явной схеме / Материалы Восемнадцатой Международной конференции по Вычислительной механике и современным прикладным программным системам (ВМСППС'2013) (Алушта, МАИ, 22-31 мая 2013 г.) С. 156-158.
7. Г.И. Баренблатт, В.М. Ентов, В.М. Рыжик Движение жидкостей и газов в природных пластах. М., Недра, 1982. - 208 с.
© А. И. Шангараева - аспирантка ИММ им. Н.И.Лобачевского К(П)ФУ, shai572@mail.ru.
© A. I Shangaraeva, post-graduate student, Institute of Mathematics and Mechanics by N. I. Lobachevsky, Kazan Federal University, shai572@mail.ru.