CC BY
95
УДК 519.63
Д. А. Мустафина, А. П. Скибин
КОНЕЧНО-ЭЛЕМЕНТНЫЙ МЕТОД КОНТРОЛЬНОГО ОБЪЕМА ДЛЯ РЕШЕНИЯ ЗАДАЧ ПОДЗЕМНОЙ ГИДРОДИНАМИКИ
Разработан вариант метода конечных элементов, позволяющий моделировать процессы фильтрации в неоднородных и анизотропных пластах со сложной геометрией. В методе используются неструктурированные сетки, содержащие треугольные элементы, что позволяет решать нелинейные задачи на грубых сетках благодаря выполнению локального закона сохранения.
В задачах моделирования процессов фильтрации рассматриваются течения, возникающие в геологических структурах сложной геометрии, таких как многопластовые коллекторы с трещинами и разрывами или призабойные зоны. Для моделирования подобных течений широко используются метод конечных элементов (МКЭ) [1] и метод контрольного объема (МКО) [2].
Основной особенностью МКО является локальное применение законов сохранения, что дает возможность прямой физической интерпретации результирующих разностных уравнений. Это привело к тому, что МКО стал предпочтительной технологией, используемой во многих коммерческих кодах (Fluent, Star-CD, CFX).
Наиболее широко распространен ячейко-центрированный вариант МКО [3]. Однако в последнее время стал применяться и вершино-центрированный вариант, что связано с применением полиэдральных контрольных объемов с существенно вершино-центрированной схемой дискретизации (STAR-CCM+) [4].
Вариант конечно-элементного метода контрольного объема (МКЭКО), разработанный Балигой [5] и описанный Пракашем [6], был успешно применен для задач гидрогазодинамики в областях со сложной геометрией [7]. Этот метод может рассматриваться как МКО с вершино-центрированной дискретизацией [8]. Многогранный контрольный объем строится вокруг узла сетки. Локальное изменение переменной внутри элемента описывается простыми кусочно-полиномиальными функциями, определенными на элементе, что позволяет получать дискретный аналог для произвольных неструктурированных сеток. Таким образом, МКЭКО комбинирует положительные черты метода контрольного объема и метода конечных элементов [9].
Форсит [10], Фунг и др. [11] использовали МКЭКО для однофазной и многофазной фильтрации, решая уравнение пьезопроводности.
Форсит применил МКЭКО для задач с локальными сеточными сгущениями, обеспечив гладкие переходы от крупной сетки к мелкой. Фунг и др. [11] внедрили МКЭКО в термальный симулятор. В работе [12] Эймард и Сониер исследовали свойства разностных схем МКЭКО. Верма и Азиз [13] разработали несколько способов создания двух- и трехмерных сеток для моделирования пластовых систем и расчетов потоков на этих сетках. Матхай и др. [14], Класс и др. [15] продемонстрировали успешное применение многосеточного метода решения матричных уравнений для многофазной многокомпонентной модели. Во многих работах [16-22] описано использование МКЭКО для двух- и трехмерных несмешивающихся течений. Особенностью данных работ является то, что поле давления определяется стандартным методом Галеркина, а поле насыщенности получается применением МКО.
В настоящей работе предлагается вариант МКЭКО на треугольных элементах (для декартовой и цилиндрической систем координат), имеющий равный порядок интерполяции скорости и давления для моделирования течения жидкости или газа в пористой среде. В предложенном методе поле давлений рассчитывается с помощью МКО, который в случае несжимаемого течения с постоянным коэффициентом проницаемости на элементе в декартовой системе координат приводит к тому же виду матриц жесткости, что и метод Петрова-Галеркина при использовании линейных функций формы [11].
Основные уравнения. Рассмотрим течение жидкости или газа в пористом двумерном теле произвольной формы с границами Г1, Г2, Г3 при следующих допущениях: жидкость является вязкой, ньютоновской и сжимаемой средой; течение описывается нелинейным законом фильтрации (уравнением Форхгеймера) [23]; для описания теплообмена используется однотемпературная модель [24].
Течение жидкости в пористом теле описывается системой уравнений, содержащей уравнение неразрывности
^ + а«(рй)=0 (1)
(где ф — пористость; р — плотность, кг/м3; Ьь — вектор скорости фильтрации, м/с; t — время, с), уравнение движения [23]
- [К] grad Р = ц [Е] ьь + р [К] [в] | Ь (2)
(здесь р — вязкость Па-е; Р — давление, Па; |ЬЬ| — модуль век-
тора скорости, м/с; [K] =
k k
™хх ™хy
kk
™ух "'yy
тензор проницаемости;
[в] =
— тензор инерционных коэффициентов; [E] — еди-
вЖЖ вЖ
_вуж вуу_
ничный тензор).
Уравнение (2) можно переписать в следующем виде:
гг = — [D] grad P,
(3)
где
[D] = [^[E] + р[К] [в]|гг| ]-1 [K] = Составляющие тензора [D] в глобальной системе координат имеют
Dyy
вид
Dra
кжу куж вуу р 1 гг | + кжж(р + kyy вуу р|гг |)
X
X
р + (к ЖЖвЖЖ + кужвжу + kxy вуж + kyy вУУ )РР|гЯ +
+ (kxy куж кжжкуу )(вжу вуж вЖЖвУУ )р |гЯ 1
1
Dжy
кжу + ( кжу куж кжж kyy ) вжу р 1 гЯ 1
X
X
р + (к жжвжж + кужвжу + кжу вуж + куу вУУ)рр|гЯ +
+ (кжу куж кжжкуу )(вжу вуж вжжвуу )р |гЯ 1
1
Dyж
куж + ( кжу куж кжж куу ) вуж р 1 гЯ 1
X
X
р + (к жжвжж + кужвжу + кжу вуж + куу вУУ)рр|гЯ +
+ (кжу куж кжжкуу )(вжу вуж вжжвУУ )р |гЯ 1
1
Dyy
кжукужАжжр|г 1 + куу (р + кжжвжжр|гЯ |)
X
X
р + (к жжвжж + кужвжу + кжу вуж + куу вУУ)рр|гЯ
+ (кжу куж кжжкуу )(вжу вуж вжжвуу )р |гЯ 1
1
Если главные оси анизотропии совпадают с осями глобальной системы координат (ж, у), тензор [Д] является диагональным:
[D] =
Du 0 0 Dv
Du =
Dv =
А + рДжж |гг 1 ^ + рвУУ |гг 1
(4)
к
к
УУ
В случае течения, подчиняющегося линейному закону фильтрации [23], составляющие тензора [в] равны нулю.
1
1
Уравнение энергии, описывающее теплообмен в пористой среде [24], имеет вид
дТ
(pcp)e— + pcpw grad T = div([Ae] grad T) + St, (5)
где cp — удельная теплоемкость, Дж/(кг-К); [Ae] =
Aexx Aexy Aeyx Aeyy
— тен-
зор эффективной теплопроводности пористой среды; Бт — мощность внутреннего источника теплоты, Дж/(м3-е); (рср)е — объемная эффективная теплоемкость, Дж/(м3-К).
Запишем обобщенное конвективно-диффузионное уравнение переноса для скалярной переменной Ф: д
— (фрФ) + div(ргo Ф) = div(ГgradФ) + БФ, (6)
где Г — коэффициент диффузии.
Заданы следующие граничные условия.
Для уравнения движения и неразрывности:
1) Р|Г1 = Рзад (постоянство давления, Па);
2) (рт,и) |Г2 = 0 (непроницаемость границы);
3) Я = (рьь,п)1Г (расход жидкости, кг/(м2с)).
Для уравнения энергии:
1) Т|г = Тзад (постоянство температуры, К, — граничное условие первого рода);
2) Чт = — ([Ае] grad Т,п)1Г2 (плотность теплового потока, Вт/м2, — граничное условие второго рода);
3) —([Аe]gradТ,п)|Гз = а(Т — Tf) (конвективный теплообмен на границе), а — коэффициент теплоотдачи, Вт/(м2К); Tf — температура окружающей среды, К.
Для уравнения переноса граничные условия аналогичны граничным условиям для уравнения энергии.
В качестве начальных условий зададим начальные распределения всех зависимых переменных: Р|4=0 = Рнач, Т|4=0 = Тнач, Ф^=0 = Фнач.
Разбиение расчетной области на контрольные объемы. Формирование контрольных объемов около узлов, находящихся в вершинах конечных элементов проводится следующим образом [5] (рис. 1): в конечном элементе определяется положение центра масс, соответствующее точке 2, затем определяются середины сторон элемента (точки 1, 3, 4), которые соединяются отрезками медиан с центром масс. Часть области, ограниченная ломаной 1-1-2-3-1, является вкладом 1-го элемента \-j-k в контрольный объем ¿-го узла; часть области ]-4-2-1-является вкладом 1-го элемента в контрольный объем j -го узла и часть к-4-2-3-к является вкладом 1-го элемента в контрольный объем ^го узла. Объединяем вклады от каждого элемента, содержащего данный
Рис. 1. Разбиение расчетной области
узел, в одну область, которая будет называться контрольным объемом для данного узла. Контрольный объем для узла г выделен штриховкой на рис. 1. Полученные таким образом контрольные объемы обладают следующими отличителными свойствами: они не перекрываются, полностью заполняя расчетную область; их границы не включают межэлементных границ; они могут быть построены на любой треугольной сетке, включающей и тупоугольные треугольники.
Интерполяция на элементе. На элементе ¡-¡-к могут применяться два вида интерполяции — линейная и ступенчатая. Линейная интерполяция: Ф = N Ф; + N Ф, + NФк, где N = N N N] — функции формы, которые для треугольного линейного элемента имеют вид, представленный в [1]. Если на элементе используется ступенчатая интерполяция, то считается, что значение функции Ф постоянно по всему контрольному объему для данного узла и равно значению функции в узле.
Случай анизотропного пласта. В метод включен общий случай анизотропии пласта. Главные оси проницаемости (ж7, у7) и теплопроводности (ж77, у77) на элементе повернуты относительно глобальных осей координат (ж, у) на углы х и в соответственно. Для удобства будем использовать локальные системы координат (X7, У7), (X77, У77), оси которых совпадают с главными осями проницаемости и теплопроводности соответственно, а начало координат совпадает с началом глобальной системы координат. Тогда тензоры [в], [К] в системах координат (ж, у), (ж7, у7), (X7,У7) принимают вид
[в ] =
вжж вуж
вжу вУУ
'вж' 0
0
вУ'
"вх' 0
0
ву'
k k xy kx' 0 kx' 0
k k kyy 0 ky'_ 0 kY'
[K ] =
В локальной системе координат (X', Y') тензор [D] запишем как
1 1
[D] =
Du 0
0
Dv
Du =
+ pßx' | w |
kx'
Dv =
+ pßy |w |
Aexx Aexy Aex'' 0 Aex'' 0
Aeyx Aeyy 0 Aey'' 0 A ey
Представим тензор эффективной теплопроводности [Ae] в системах координат (x,y), (x'',y''), (X'', Y''):
[Ae] =
Равный порядок интерполяции скорости и давления. Трудности метода равного порядка интерполяции связаны с тем, что если давление интерполируется с использованием линейной функции формы, то только разность давлений между чередующимися (каждым вторым) узлами включается в полную систему уравнений. Следовательно, разностные уравнения не могут обнаружить разницу между равномерным и шахматными полями давлений. В конечно-разностных формулировках МКО предотвращение шахматного поля давлений осуществляется за счет использования смещенных сеток [2].
В конечно-элементных методах шахматное поле давления исключается посредством применения разных сеток конечных элементов для определения скорости и давления [25]; получения давления с использованием уравнения Пуассона, как это сделано Шнейдером и др. [26]; использования некоторых специальных методик для фильтрации нереальных колебаний давления [27]; использования формулировок через штрафные функции [28].
В данной работе рассматривается модификация метода равного порядка интерполяции скорости и давления [6], позволяющего избежать шахматного поля давления, для уравнения Форхгеймера. Основная идея метода [6] заключается в том, что массовый расход, рассчитанный по скоростям У, определяемым в узлах, заменяется расходом, рассчитанным по скоростям У, ассоциированным с 1-м конечным элементом. Эта идея схожа с подходом использования смещенных сеток. Поле скорости У = йвх> + йеу определяется на каждом элементе и имеет разрывы в его узлах (ех>, £у> — единичные орты системы координат (X' ,У').
Компоненты поля скорости У определяются следующим образом: дР
U = - Du
X'
' dN' {P}; v = v dP "dN"
= —Du D д
i dX'_ dY' = -Dv i dY'_
{P },
dP
где
dX'
dP
и ^ттт — постоянные значения градиентов давления на /-м
конечном элементе г——Значения коэффициентов и определяются в соответствии с выражениями (7).
Вывод дискретного аналога для уравнения неразрывности.
При решении уравнения неразрывности сделаны следующие допущения: направление главных осей анизотропии на элементе постоянно; плотность, теплопроводность, теплоемкость, пористость непрерывны на элементе и изменяются линейно. Запишем уравнение неразрывности (1) в виде
фдр + (рйи) = 0.
Применим локальный закон сохранения массы к каждому контрольному объему. Тогда для расчетной области, состоящей из М контрольных объемов, имеем следующую систему уравнений:
/ ^+у
Ух Ух
I +/
div (^piuj dV = 0; div (pwU ) dV = 0;
V2
J Ф^У + ^ div (рйТ) ¿V = 0.
^ Ум Ум
где У1, У2,... Ум — контрольные объемы.
Эта система после вычисления интегралов переходит в систему разностных уравнений. Интегрирование будем проводить по вкладам элементов в контрольный объем. Так, для контрольного объема г-го узла получим
L /»о L л
, ^+£
14 V,
Е
div (piu) dV = 0.
(8)
где Уц — вклад от /-го элемента в контрольный объем г-го узла; Ь — число элементов, окружающих г-й узел.
Используя теорему Гаусса-Остроградского, находим:
/ div(pгй)dУ — / (рй,п)^£ = / (Р,
где . = рй = — р;[А] grad Р — массовый поток; п — вектор внешней нормали. Далее индекс / у величины р и [Д] опустим.
Уравнение (8) содержит нестационарный и дивергентный члены. Суммируя их, получим дискретный аналог уравнения. Рассмотрим нестационарный член уравнения (8)
Г = *,„^ ,
Уг1
где индекс "*" соответствует значению на предыдущем шаге по времени.
В этом случае для ф и р применена ступенчатая интерполяция. В случае слабосжимаемой жидкости получим
[ др [ дР Рг — Р*
= ФйРхс№ = ФйРаСт 'Аг г АУй,
Vil
Vil
где cf =
1 dp PdP
— сжимаемость жидкости, 1/Па.
т
Рассмотрим дивергентный член. Интегрирование по внешней поверхности следует проводить только для граничных конечных элементов, так как для внутренних элементов интегралы взаимно уничтожаются. Определим интеграл по внутренней поверхности для узла i:
J(F,n)dS = J(F,n)dS + J(F,n)dS.
Sil 12 23
Вследствие инвариантности скалярного произведения имеем - dP dP
(F,n)(x,y) = (F,n)(x',y') = (F,n)(X',Y') = PD'UdXnX' — PD'V ~gY'',
где Du, Dv определяются соотношением (7).
Для интерполяции давления на элементе используется линейная интерполяция:
J(F,n)dS = J (F,n)dS = J(F,n)dS + j{F,n)dS,
Sil i—l—2—3—i 12 23
(F,n)dS = — ^nx' f pDu dS — dPluy' [ pDv dS—
dY'
Sii
12
12
— Шux' [ pdu ds — dpUY' f Pdv ds =
= (pDu)12(Y/ — Y2')
dY'
23
' dN
dx'
23
{P} + (pDv )12(X2 — X1)
dN
W'
{p }+
+ (pDu) (Y' — Уз')
dN
dx7
{p } + (pDv )23(X3 — X7)
dN
dY7
{p },
где
(pDu)12 = -1 f pDu dS; (pDv)12 = 1 f pDv dS; П2 J 112 J
12 12
/12 = V(X1 — X2 )2 + (У1' — Y2')2;
(pDu)23 = -1 f pDu dS; (pDv)23 = 1 f pDv dS; '23 J '23 J
23
23
/23 = >/(X2 — X3 )2 + (Y7 — Y3')2.
Затем, проведя аналогичные вычисления для вкладов от других элементов, получим следующий дискретный аналог уравнения неразрывности для г-го узла
(£ аПь) Pp = £ an>Pnb + dP.
Алгебраические аспекты вывода этого уравнения, аналогичные построению дискретного аналога для уравнения теплопроводности, изложены в работе [29].
Определение поля скорости в узлах сетки рассмотрим на примере узла г. Запишем в разностной форме уравнения движения на элементе:
dP
«X, = dX
dP
; ^^ = —-d ^
й1 = y^
« w + ;
го = их' ех' + -Уу/еу/.
г1 г1
Тогда компоненты скоростей на элементе в глобальной системе координат следующие:
иц = (ггц, еХ); = (гогЬ еу).
После интегрирования по объему получим значения компонент скорости в узле:
/ «ii dV
7_1
■Vi,
1
«i =
£ AV^u^;
£ i Vii dV
7_1
Vi =
-Vi,
1
Vi " ' ' Vi Vi
i i i=i i i i=i
£ AViiVii.
где V = ^^ A Vi — контрольный объем i-го узла.
i
Особенности дискретного аналога уравнения переноса. Уравнение переноса (6) содержит нестационарный, конвективный, диффузионный, а также источниковый члены [2]. Рассмотрим только особенности дискретизации конвективного члена, так как процедура вычисления дискретных аналогов остальных членов описана выше и приведена в работе [29].
Так, для вклада 1-го элемента в контрольный объем ¿-го узла имеем
[ (Р • Ф,и^Б = ( [Р • + [ (Р • Ф,п^Б,
Sii
12
23
где
F • Ф,nj dS = Ф12у (F,n)dS; 12 12
F • Ф, nj dS = Ф23 у (F, nj dS.
23 23
Для функции Ф используется противопоточная интерполяция [30]:
Ф12 = <
Ф23 =
Ф, Ф,
Ф»,
Фк,
(F,n)dS > 0 : г ^ j;
12
(F,n)dS < 0 : j ^ г,
12
(F,n)dS > 0 : г ^ к;
23
(F,n)dS < 0 : к ^ г.
23
Применение метода. Для иллюстрации эффективности разработанного метода решен ряд тестовых задач: изотермическая стационарная фильтрация жидкости; изотермическая стационарная фильтрация газа с нелинейным законом фильтрации через пористый цилиндр; задача пробоотбора.
Изотермическая стационарная фильтрация жидкости. Рассмотрим стационарное течение жидкости плотностью р = 1000 кг/м3, вязкостью р = 0,01 Па-с через изотропный пласт радиусом ЯЕ = 3,1 м с проницаемостью к = 10-13 м2 и пористостью ф = 0,2 к скважине радиусом Ящ = 0,1 м. Жидкость фильтруется под действием перепада давлений: пластовое давление Рщ = 25 МПа, давление на забое
РЕ = 30 МПа. Необходимо найти распределения давления и скорости жидкости. Задача решается в цилиндрической системе координат. Аналитические выражения для давления и скорости следующие:
P=
ре- pw
w
ln (RE/RW)
ln r +
Pw ln Re - Ре ln R
w
kdP
v = —-
ln (Re/Rw) 1 k Ре — Pw
р дг тр 1п^е/RW) Численное решение получено с использованием равномерной и неравномерной (сгущенной к скважине) сеток 11x11. Сравнение результатов показало хорошее соответствие между аналитическим и численным решениями (табл. 1).
Изотермическая стационарная фильтрация газа с нелинейным законом фильтрации через пористый цилиндр. Рассмотрим течение воздуха вязкостью р = 18,040-6 Па-с через пористый цилиндр (внутренний радиус Т1 = 0,03 м, внешний радиус т2 = 0,04 м) с проницаемостью к = 1Д5432М0-11 м2 и пористостью ф = 0,35. Давление на внутренней поверхности Р1 = 1 МПа, на внешней — Р2 = 0,1 МПа. Температура воздуха постоянна и равна Т = 273 К. Воздух считается идеальным газом. Инерционный коэффициент в = 6,308404 м-1. Необходимо определить распределение давления по радиусу и расход воздуха.
Таблица 1
Давление P, МПа (равномерная сетка 11x11) Давление P, МПа (неравномерная сетка 11x11)
r Точное Численное Погреш- r Точное Численное Погреш-
решение решение ность, % решение решение ность, %
0,1 25,000 25,000 0,00 0,100 25,00 25,00 0,00
0,4 27,018 26,840 0,66 0,184 25,89 25,87 0,07
0,7 27,833 27,700 0,48 0,316 26,68 26,65 0,09
1,0 28,353 28,250 0,36 0,496 27,33 27,31 0,08
1,3 28,735 28,650 0,29 0,724 27,88 27,86 0,08
1,6 29,037 28,980 0,20 1,00 28,35 28,33 0,08
1,9 29,287 29,240 0,16 1,324 28,76 28,75 0,04
2,2 29,501 29,470 0,10 1,696 29,12 29,11 0,04
2,5 29,687 29,670 0,06 2,116 29,44 29,44 0,01
2,8 29,852 29,840 0,04 2,584 29,73 29,73 0,02
3,1 30,000 30,000 0,00 3,100 30,00 30,00 0,00
Расход G, кг/(м 2-с) Расход G, кг/(м 2-с)
9,144E-02 9,728E-02 6,38 9,144E-02 9,276E-02 1,44
Таблица 2
V — Vi Г2 — Vi Давление, МПа (численное решение) Давление, МПа Погреш-
Сетка (точное решение) ность, %
11 х 11 21 х 21 41 х 41 61 х 61
1,0 1,000 1,000 1,000 1,000 1,000 0
0,9 3,207 3,137 3,023 3,008 2,976 7,76-1,08
0,8 4,373 4,237 4,169 4,156 4,133 5,81-0,56
0,7 5,265 5,137 5,094 5,085 5,069 3,87-0,32
0,6 6,045 5,938 5,910 5,905 5,893 2,58-0,20
0,5 6,766 6,679 6,661 6,658 6,649 1,76-0,14
0,4 7,449 7,380 7,369 7,368 7,361 1,20-0,10
0,3 8,107 8,054 8,049 8,049 8,044 0,78-0,06
0,2 8,749 8,711 8,710 8,710 8,707 0,48-0,03
0,1 9,380 9,358 9,359 9,359 9,357 0,25-0,02
0,0 10,000 10,000 10,000 10,000 10,000 0
Расход воздуха G, кг/(м2-с)
4,592 4,785 4,776 4,735 4,703 2,36.. .0,68
Имеем выражения для давления
,2 1 »„„гп^г ß^^г — Г!
Р = Р2 = - увЯТ 1п--- С'ЯТ-
2 к Т1 8 тт1
и расхода воздуха
в = 2^Т2Т 1п(т2/Т1) + I СР.Т2Т1 1п(т2/Т1Л 2 + Т2Т1(Р2 - Р|)
= вк(Т2 - Т1) + V V 2вк(Т2 - Т1) ) + 2вЯТ(Т2 - Т1)'
Численное решение проводим на сетках с числом узлов 11x11, 21x21, 41x41, 61x61 расчетной области. Результаты хорошо согласуются с аналитическим решением (табл. 2).
Задача пробоотбора. После процесса бурения скважины радиусом Яш = 0,1 м пласт радиусом Яе = 3,1 м и высотой Н = 10 м заполнен буровым раствором на глубину На = 0,1 м. Пласт имеет пористость ф = 0,3 и проницаемость к = 0,5^ 10-12 м2. Нефть отбирается через поверхность отбора высотой Н,и, = 0,01125 м, куда попадает вместе с буровым раствором под действием перепада давлений (пластовое давление Ре = 20 МПа, давление на забое Р,и, = 15 МПа). Свойства нефти и бурового раствора считаются одинаковыми: вязкость
Рнефть = Рбур.раств = Р = 0,01 Пас Рнефть = Рбур.раств = Р = 1000 кг/м3.
Фазовые проницаемости предполагаются линейными функциями насыщенности кнефть = £нефть, кбур.раств = ¿бур.раств. Необходимо определить расход нефти через поверхность отбора и среднеинтегральную
Рис. 2. Схема расчетной области
концентрацию нефти на поверхности отбора в зависимости от времени.
Схема расчетной области приведена на рис. 2. Математическая модель задачи содержит уравнения неразрывности для обеих фаз
div
k | РнеФть кнефть + рбур.раств кбур.раств | grad P \ Мнефть Мбур.раств J
= о
(так как $нефть + ^бур.раств 1? кнефть ^нефть? кбур.раств ^бур.раств и
свойства нефти и бурового раствора считаются одинаковыми? предыдущее уравнение неразрывности может быть переписано в виде gгad Р] = 0) и уравнение для концентрации нефти
д
дЦ (ФРСнефть) + ¿IV(Р^ЯСнефть) = 0.
Так как плотности постоянны? то массовые концентрации равны на-сыщенностям.
Начальные и граничные условия для уравнения неразрывности имеют вид:
Р (г, г, 0) = Ре, г е (0,Н), г е (Яы, Яе); Р (Яы, г, ¿) = Pw, г е (0, Hw); Р (Ве,г,Ь) = Ре, г е (0,Н); Р (г,Н,г) = Ре, г е (Вы,Ве);
=0, Г е (Bw.Be); дР<Я^ = 0, г е (Н„Н).
Начальные и граничные условия для уравнения концентрации нефти:
Снефть (г, 0) = 1, г е (Яш + На,ЯЕ), г е (0, Н);
Снефть (г, 0) = 0, г е (Яш,ЯЖ + На), г е (0,Н);
Снефть (г, Н, ¿) = 0, г е (Яш, Яш + На);
Снефть (г,Н,£) = 1, г е (Яш + На,Яе);
Снефть (Яе,г,£) = 1, г е (0, Н);
дС"*д^М) =0, г е (Яш,Яе); дСнефть =0, г е (0,Нш).
дг
Для демонстрирации эффективности предлагаемого метода задачу решали на двух сетках. Сравнение результатов проводили с коммерческим программным обеспечением ЕсНр8е-100. Сетка 1 имеет прямоугольную топологию и 2601 узел (рис.3,а), она же используется в ЕсПр8е-100. Сетка 2 имеет сферическую топологию и 1927 узлов (рис. 3, б). Сгущение к поверхности отбора для обеих сеток приведено на рис. 4.
Среднеинтегральная насыщенность нефти на поверхности отбора в зависимости от времени приведена на рис. 5. Результаты численного расчета хорошо согласуются с результатами, полученными с помощью ЕсНр8е-100. Значения расхода нефти через поверхность отбора следующие: 7,78-10-2кг/(м2•с) (сетка 1); 7,69-10-2кг/(м2•с) (сетка 2); 7,16-10-2кг/(м2•с) (ЕсИр8е-100).
Выводы. Разработан вариант МКЭКО с равным порядком интерполяции скорости и давления для решения задач однофазной фильтрации. Оценка метода проводилась с помощью ряда тестовых задач:
Рис. 5. Зависимость насыщенности нефти на поверхности отбора от времени:
сплошная и штриховая кривые — расчет по МКЭКО (сетки 1 и 2 соответственно); о — расчет с использованием ПО "ЕсНр8е-100"
изотермическая стационарная фильтрация жидкости, изотермическая стационарная фильтрация газа с нелинейным законом фильтрации через пористый цилиндр, задача пробоотбора. Получено хорошее соответствие результатов численного моделирования с аналитическими решениями и решениями, полученными с помощью коммерческого программного обеспечения.
Авторы благодарят компанию "Шлюмберже" за предоставленную возможность проведения исследований на базе Московского научно-исследовательского центра.
СПИСОК ЛИТЕРАТУРЫ
1. Зенкевич О.,Морган К. Конечные методы и аппроксимация. - М.: Мир, 1980.
2. Патанка р С. Численные методы решения задач теплообмена и динамики жидкости. - М.: Энергоатомиздат, 1984. - 152 с.
3. Bertolazzi E., Manzini G. A cell-centered second-order accurate finite volume method for convection-diffusion problems on unstructured meshes // Mathematical Models and Methods in Applied Sciences. - 2004. - 14(8). - P. 12351260.
4. P e r i c M. Simulation of flows in complex geometries: new meshing and solution methods / NAFEMS Seminar; Simulation of Complex Flows (CFD)-Applications and Trends. Niedernhausen/Wiesbaden, Germany, 2004.
5. B a l i g a B. R., Patankar S. V. A new finite-element formulation for convection-diffusion problems // Numerical Heat Transfer. - 1980. - No. 3. - P. 393409.
6. P r a k a s h C, P a t a n k a r S. V. A control volume-based finite-element method for solving the Navier-Stokes equations using equal-order velocity-pressure interpolation // Numerical Heat Transfer. - 1985. - No. 8. P. 259-280.
7. Reyes M., R i n c o n J., D a m i a J. Simulation of turbulent flow in irregular geometries using a control-volume finite-element method // Numerical Heat Transfer.
- 2001. Part B(39). - P. 79-89.
8. T a y l o r G. A., B a i l e C., C r o s s M. A vertex-based finite volume method applied to non-linear material problems in computational solid mechanics // International Journal for Numerical Methods in Engineering. - 2003. - No. 56. -P. 507-529.
9. M c B r i d e D., C r o f t T. N., C r o s s M. A coupled finite volume method for the solution of flow processes on complex geometries // lnt. J. Numer. Meth. Fluids. -2007. No. 53. - P. 81-104.
10. Forsyth P. A. A control-volume, finite-element method for local mesh refinement in thermal reservoir simulation // SPERE, (Nov. 1990). - P. 561-566.
11. F u n g L. S. K., H i e b e r t A. D., N g h l e m L. X. Reservoir simulation with a control-volume finite element method // SPERE (Aug. 1992). - P. 349-357.
12. EymardR., SonierF. Mathematical and Numerical properties of control-volume, finite-element scheme for reservoir simulation // SPERE (Nov. 1994). -P. 283-289.
13. Verm a S., Azi z K. A control volume scheme for flexible grids in reservoir simulation // Paper SPE 37999 presented at Reservoir Simulation Symposium. -Dallas, June 8-11, 1997.
14. M a t t h a i S. K., M e z e n t s e v A., B e l a y n e h M. Control-volume finite-element two-phase flow experiments with fractured rock represented by unstructured 3D hybrid meshes // Paper SPE 93341 presented at Reservoir Simulation Symposium.
- Houston, January 31-February 2, 2005.
15. C l a s s H., H e l m i g R., B a s t i a n P. Numerical simulation of non-isothermal multiphase multicomponent processes in porous media. 1 An efficient solution technique // Advances in Water Resources. - 2002. - No. 25. - P. 533-550.
16. Hub e r R., Helmi g R. Node -centered finite volume discretizations for the numerical simulation of multiphase flow in heterogeneous porous media // Computational Geosciences. - 2000. - No. 4. - P. 141-164.
17. Durlofsky L. J. A triangle based mixed finite-element-finite volume technique for modeling two phase flow through porous media // Journal of Computational Physics. - 1993. - No. 105. - P. 252-266.
18. GottardiG., Dall'OlioD. A control-volume finite-element model for simulating oil-water reservoirs // J. Pet. Sci. Eng. - 1992. - No. 8. - P. 29-41.
19. C o r d a z z o J., M a l i s k a C. R., S i l v a A. F. C., H u r t a d o F. S.V. The element-based finite volume method applied to petroleum reservoir simulation // IBP29404 presented at Rio Oil & Gas Expo and Conference. - Rio de Janeiro, October 4-7, 2004.
20. GeigerS., Roberts S., MatthaiS., ZoppouC.,Burry A. Combining finite volume and finite element methods for efficient multiphase flow simulations in highly heterogeneous and structurally complex geologic media // Geofluids. - 2004.
- No. 4. - P. 284-299.
21. N a j i H. S., K i n g A., K a z e m i H. A fully implicit, three-dimensional, two-phase, control-volume Finite element model for the simulation of naturally fractured reservoirs // Paper SPE 36279 presented at the 7th Abu Dhabi international petroleum Exhibition and Conference. - Abu Dhabi, October 13-16, 1996.
22. F u Y., Yang Y. -K., D e o M. Three-dimensional, three-phase discrete-fracture reservoir simulator based on control volume finite element (CVFE) formulation // Paper SPE 93292 presented at SPE Reservoir Simulation Symposium. - Houston, January 31-February 2, 2005.
23. HuangH.,AyoubJ. Applicability of the Forchheimer equation for non-Darcy flow in porous media // Paper SPE 102715 presented at Annual Technical Conference and Exhibition. - San Antonio, September 24-27, 2006.
24. Ч е к а л ю к Э. Б. Термодинамика нефтяного пласта. - М.: Недра, 1965. - 240 с.
25. H u y a k o m P. S., T a y l o r C., Lee R. L., C i r e s h o P. M. A comparison of various mixed-interpolation finite elements in the velocity-pressure formulation of the Navier-Stokes equations // Comput. Fluids. - 1978. - No 5. - P. 25-35.
26. S h n e i d e r G. E., R a i t h b y G. D., Y o v a n o v i c h M. M. Finite-element solution procedures for solving the incompressible Navier-Stokes equations using equal order variable interpolation // Numer. Heat Transfer. - 1978. - No 1. - P. 433451.
27. S a n i R. L., G r e s h o P.M., Lee R. L., G r i f f i t h s D. F. The cause and cure of the spurious pressures generated by certain FEM solutions of the incompressible Navier-Stokes equations. Part I // Int. J. Numer. Meth. Fluids. - 1981. - No 1. -P. 17-30 (Also see Part II. - P. 171).
28. Heinric h C., Marshall R. S. Viscous incompressible flow by a penalty function finite element method // Comput. Fluids. - 1981. - No. 1. - P. 73-83.
29. С к и б и н А. П., Ч е р в я к о в В.В.,Югов В. П. Метод конечных элементов, основанный на интегрировании по контрольному объему для двумерных нестационарных эллиптических задач // Изв. РАН. Сер. "Энергетика". - 1995. -№ 1.
30. A z i z K., S e 11 a r i A. Petroleum reservoir simulation, London: Applied Science Publishers. - 1978.
Статья поступила в редакцию 10.10.2007
Александр Петрович Скибин родился в 1963 г., окончил МВТУ им. Н.Э. Баумана в 1986 г. и МГУ им. М.В. Ломоносова в 1988 г. Канд. техн. наук, доцент кафедры "Теплофизика" МГТУ им. Н.Э. Баумана, старший научный сотрудник научного центра компании "Шлюмберже". Автор более 70 научных работ в области разработки и применения методов математического моделирования для инженерных приложений.
A.P. Skibin (b. 1963) graduated from the Bauman Moscow Higher Technical School in 1986 and Lomonosov Moscow State University in 1988. Ph. D. (Eng.), assoc. professor of "Thermal Physics" department of the Bauman Moscow State Technical University, senior worker of scientific center of the company Schlumberger. Author of more than 70 publications in the field of development and use of methods of mathematical simulation for engineering applications.
Дарья Александровна Мустафина родилась в 1985 г., студентка кафедры "Теплофизика" МГТУ им. Н.Э. Баумана. Специализируется в области методов математического моделирования процессов тепло- и массообмена.
D.A. Mustafina (b. 1985) — student of "Thermal Physics" department of the Bauman Moscow State Technical University. Specializes in the field of methods of mathematical simulation of processes of heat and mass exchange.
