Содержание к диссертации
Введение
1 Введение 7
1.1 Описание явления 7
1.2 Актуальносіь рассматриваемой задачи 7
1 3 Объект и предмет исследования 9
1 4 Состояние проблемы и методы описания 9
14.1 Исследование фроніа в однородных средах 9
1 4,2 Исследование фроніа в неоднородных средах 12
14 3 Стохастический анализ двухфазных потоков 13
1.4.4 Численное моделирование двухфазных потоков 21
1.4.5 Методы ренормализационной группы и апскейлинга . 23
1.5 Цель и задачи работы 25
1.5.1 Стохастическое моделирование 25
1.5.2 Численное моделирование 26
1.5 3 Моделирование устойчивого фронта 26
1 б Методы исследований 28
1.7 Научная новизна 29
1 8 Основные положения, выносимые на защиту 30
19 Апробация работы 30
1.10 Основные публикации по теме диссертации 31
1.11 Краткий обзор содержания работы 32
2 Стохастическое моделирование распространения фронта 35
2.1 Определяющие уравнения 35
2.2 Уравнение для формы фронта 39
2.3 Скороеіь фильтрационного потока 44
2 4 Уравнение для флуктуации формы фронта 50
3 Исследование статистических характеристик фронта 59
3.1 Корреляционная функция и дисперсия флуктуации формы фронта 59
3.2 Дисперсия продольных скоростей на фронте вьпеснения 63
3.3 Вариограмма продольных смещений фронта вытеснения 65
3 4 Средняя насыщенность и дисперсия насыщенности 70
3.5 Сравнение с результатами вычислиіельиьіх экспериментов 80
4 Численный метод быстрого моделирования распространения фронта 93
4.1 Основная идея меюда 93
4.2 Общее уравнение 95
4 3 Численная схема 97
4 3.1 Дискретизация однофазной неоднородной задачи 97
4 3.2 Дискретизация члена, учитывающего неоднородное в среды 100
4 33 Дискретизация в просгранстве Фурье-образов 100
4 3 4 Дискретизация по времени 102
4.4 Результаты и сравнение 105
5 Заключение 112
5.1 Обюр полученных результатов 112
5.2 Возможные направления дальнейших исследований 114
А Приложения 116
АЛ Вычисление насьіщенносіи на задней сгороие поверхности фронта 116
А.2 Использование функции Грина 117
A 3 Вычисление ишеграла от спектральной плогносш для корреляционной функции 120
А.4 Вычисление интеграла от экспоненциальной функции 123
А.5 Вычисление интеграла ог спектральной илоіііосіи для вариограммы 124
В Основные обозначения 125
Список литературы
- Состояние проблемы и методы описания
- Уравнение для формы фронта
- Дисперсия продольных скоростей на фронте вьпеснения
- Дискретизация члена, учитывающего неоднородное в среды
Введение к работе
Работа посвящена моделированию несмешивающихся двухфазных течений в случайно-неоднородных пористых средах в применении к задаче вытеснения. Эта задача актуальна в нефтяной индустрии при использовании путей повышения неф геоидами.
Состояние проблемы и методы описания
Последующие исследования задачи вьпеснеиия представляли собой анализ сіойчивосіи плоского фроша по отношению к малым возмущениям в двумерной однородной пористой среде. В частности, к ним относяіся классические работы Сафмена и Тейлора [5] и Чуока с соавюрами [б]. Авторы рас сматривали процесс, при котором одна жидкосіь полное і ью вытесняла дру-іую не смешиваясь с ней при эюм толщина переходной зоны возле границы раздела как правило мала по сравнению с макроскопическими масштабами сисіемьі. Исследуя границу раздела между жидкостями одинаковой плотности и используя условие непрерывности нормальных компонент скорости и давления на фронте, авторы пришли к следующему заключению- граница раздела устойчива по о і ношению к малым возмущениям, если более вязкая жидкосіь вытесняет менее вязкую, и наоборот — граница раздела неустойчива к малым возмущениям, если менее вязкая жидкость вытесняет более вязкую. Иными словами, вводя параметр т. представляющий собой отношение вязкости вытесняемой жидкое і и к вязкости вытесняющей, граница раздела усюйчива по отношению к малым возмущениям, если т 1 и неустойчива, если т 1. Если фронт неустойчив, то малые возмущения растут со временем, образуя сложные іеомеїрические сіруктурьі. коюрые получили название «языков» («fingers» в западной литературе) 7]. Для экспериментальною исследования таких структур используеіся усіановка. называемая ячейкой Хеле-Шоу. Эта усіановка представляет собой две прозрачные ило-скоиараллельные пластины с малым зазором 6 между ними, заполненным жидкостями. Такая система позволяет моделировать однородную пористую среду, проницаемосіь которой полагается равной 62/12. и непосредственно наблюдать течения.
Поскольку в рассматриваемой Сафменом. Тейлором и Чуоком с соавторами модели предполагалось чю одна жидкосіь вытесняет другую без остатка, то не принимались во внимание эффекты относительной проницаемости, имеющие месю в реальных пористых средах [8]. Суть отих эффектов заключается в том. что проницаемость данной среды для данной жидкости зависит от насьіщенносіи. На практике, кривые относительной проницаемости определяю іся экспериментально и эта зависимость часто оказывается нелинейной. Подобные эффекты приводят к тому, что при решении методом характеристик описывающею эволюцию насыщенности гиперболического уравнения оказывается, что решение имеет разрывы, на которые ранее обратили внимание Баклей и Левереті. Такое обобщенное решение получило название решения Баклея-Леверетіа [4]. Таким образом, в процессе вытеснения можно выделиіь область постепенного уменьшения насыщенности (также называемую «зоной разрежения») и скачок насыщенности, с которым в дальнейшем и будем ассоциировать фронг.
Необходимоеіь учеы влияния зависимое:и характеристик двухфазною течения от насьіщенносіи на устойчивость фронта вытеснения была указана в работе Хагорта [9]. В этой работе автор провел анализ устойчивости фроніа вытеснения для модели Баклея-Леверетта принимая во внимание эффекты относительной проницаемости и показал, что в этом случае устойчивость фроша по отношению к малым возмущением определяется фронтальным отношением иодвижносіей М/ = A(5i)/A(S,2) (определение подвижности см. в разделе 2.1. формула (2.1.12)). где S\ — значение насыщенносш на задней сюронс поверхности фронта а.% — значение насыщенности на передней стороне поверхности фронта Если М/ 1. то граница раздела усюйчива по 01 ношению к малым возмущениям. Если Mj 1. то малые возмущения границы раздела будут расти со временем, приводя к образованию языков. При использовании квадратичных функций для зависимости относительной проницаемости от насыщенности для моделирования двухфазных потоков (см. раздел 2.1. формула (2.1.13)) это приводит к тому, что критерий устойчивости для вязкости изменяется: фронт устойчив по отношению к малым возмущениям, если т 3 и неустойчив, если т 3 (в отличии от т 1 и т 1 в исследованиях Сафмена-Тейлора и Чуока с соавторами). Подробное исследование устойчивости течений с различными ігачальньїми распределениями насыщенности проведено в работе [10]. но для интересующих нас случаев оно сводится к условиям, полученным в работах Сафмена-Тейлора, Чуока с соавторами и Хагорта. Таким образом, учет влияния насыщенности на картину течения оказался весьма существенным.
Отметим гакже анализ устойчивости, проведенный в работе Кинга и Ду-наевскою [11]. в которой в качестве невозмущенної о сосюяния использовалось решение Баклея-Леверегта и было получено уравнение, описывающее временную эволюцию флуктуации фронта вытеснения. Полученный ими критерий устойчивости совпал е результатами Хагорта.
Уравнение для формы фронта
Предметом исследования будет задача вытеснения в случайно-неоднородной среде в пространстве произвольной размерности d (в реальности d — 1,2,3. но случай d — 1 являє і ся тривиальным).
В дальнейшем введем следующие обозначения. Совокупность пространственных координат г (с/-мерный вектор) имеет компоненты {xt} (латинские индексы принимают значения і = 1...(/). где ось х = Xi направлена вдоль усредненного направления распространения фронта: компонент ор-j 01 опальної о к оси Х\ (d — 1)-мерного вектора у обозначим через ха. где греческие индексы принимают значения (a = 2...(1). В этих обозначениях вектор ірадиеша имеет компоненты {дг} = {ді,да} а во-шикающий при выполнении преобразований Фурье волновой векюр Q = {#i,q} имеет компоненты
Пусть уравнение поверхносіи разрыва (поверхности фронта вытеснения) имеет вид у {т,у, t) = 0; Это уравнение может быть разрешено относительно х и переписано в виде ip{x,y,t) = x-h{y,t) = 0.
Функция h(y,t) описывает форму фронта Значению р 0 соотвеїствуег область, занимаемая первой жидкосіью (водой), а р 0 — область, занимаемая второй жидкосіью (нефтью). Передней и задней сторонам граничной поверхносіи соответствуют р = +0 и ір = —0.
Рассмотрим цилиндрический обьем U пористой среды с площадью поперечною сечения Е. Пусть скачок насыщенности проходит чере І эют обьем по нормали к Е так. чю поверхность скачка параллельна основаниям цилиндра и рассюяния от оснований цилиндра до скачка равны /. Условие сохранения массы вытесняющей жидкости в элементе имеет вид dtj 4 SdQ. + / uwndE = 0, где иип — нормальная к границе раздела компонент скоросіи вытесняю1 щей жидкое і и. Положим, что S\ и 6 2 — значения насыщенностей на задней и передней сюронах поверхности раздела (при моделировании двухфазных течений с учеюм эффектов оіносиїельньїх проницаемостей S\ представляет собой значение насыщенности на задней сюроне поверхности фронта, a , — значение насьіщешюсіи изначально присутствовавшее в пористой среде). Тогда изменение обьема вытесняющей жидкое і и со временем можно вырази хь следующим сооїношеїшем (см. [79]) jt J })SdU = t {Sv - S2)HVn + 0(Z/R2), где Vn — скорость перемещения скачка no нормали к нему. R — радиус кривизны поверхности скачка. Разность иоюков вытесняющей жидкости через сечения, параллельные поверхносш скачка, равна [uwn\ — ишпг)Е. где индексами 1 и 2 обозначены нормальные компоненты скоросіи вытесняющей жидкости па задней и передней сторонах поверхносіи раздела соотвеїсгвен-но Касательные сосіавляющие скоросіи исчезающе малы при / - 0. Тогда условие сохранения жидкое і и при стягивании цилиндра обьема ft к участку поверхносіи примет вид Vn да - s2)
Нормальная сосіавляюіцая суммарной скоросіи фильграции при переходе через поверхность разрыва должна быгь непрерывна, следовательно Uon + Uwn = Un\ = Un2 = Un, где ип — нормальная к поверхносіи сосіавляющая суммарной скорости фильтрации. Тогда, используя (2.1.17). получаем v _unF(S1)-F(S2) ф b\ — 2 С ді)угой сюроны. при движении границы раздела имеем at где V — век юр скорости границы раздела Это уравнение можно переписать в следующем виде v = -dMv pyl = - т р I2-2-2) а нормальная составляющая скорости определяв і ся как Vn = V . n, (2.2 3) іде n = - = l ={1, dah{y, t)}. (2.2.4) VH yi + I My,OI2 Тогда, подставляя (2.2.4) и (2.2.2) в (2.2 3). получаем Vn = - L = №1 (2.2.5)
Окончательно учитывая что un = u n и подсіавляя (2.2 4) и (2.2.5) в (2.2.1) приходим к уравнению lf( 1 FjSj) - F(S2), _ ,, Ф Ьі — i 2 = — {щ - uadQh)\v=o, (2.2.G) Щ где со = -j 5 с = 7F Wi)- i2-2-7)
Здесь «о — направленная вдоль оси х средняя скорость. Здесь и далее по тексту принимается правило суммирования по повюряющимся индексам.
Нас интересует уравнение, описывающее флуктуации фроша вытеснения, которое следует из (2.2.6). В линейном приближении необходимо сохранить член нулевого и первою порядка по h, те. предеіавить уравнение (2.2 6) в виде dth = C-h + f, (2.2.8) где С — некоторый линейный дифференциальный оператор и / — некоторая случайная функция.
Из (2.2 6) видно что для нахождения эгого уравнения в формуле (иДу ) = («і — uadah) следует оставить только член нулевою и первого порядка но h для щ и нулевою порядка в выражении для иа.
Если с помощью уравнений (2.1.16) и (2.1.9) найти скорость фильтрации и как функцию 5 и иодставиїь в (2.2 6). то получится уравнение для функции h(y,t). задающей форму граничной поверхности раздела. Принимается, чю значения нлсыщенностей но обе стороны границы раздела 5i и 5г являю і -ся заданными консіантами. Задача вычисления скорости фильтрационної о потока рассматривается в следующем параграфе.
Дисперсия продольных скоростей на фронте вьпеснения
Приведем выражение, полученное для средней насыщенное і и и дисперсии насыщенности, к форме, удобной для вычислений и графических построений. Для оюю сосіавим безразмерную величину ,/, называемую безразмерным временем, следующим образом U = t , (3.4.6) где L — характерный макроскопический масштаб системы. Если речь идет об ограниченной обласіи. то L — ее длина в направлении х. Предполаїается. чю пористость ф = 0,4. В трехмерном случае была получена связь между дисперсией смещений фронта Во и парамеїрами экспоненциальной корреляционной функции логарифмов проницаемости К = Коехр(-г/а) Во = l-I(3,A)KQa2 = 1/(3, A)K L2 После ал і ебраических манипуляций, используя (2.2.7). получаем 2 2 \ук„1(3,А))
В последнее выражение входят еще два безразмерных параметра. Первый из них x/L харакіеризуег оіносиїельное положение точки вдоль оси х. Второй парамеїр a/L представляет собой отношение длины корреляции поля логарифмов проницаемости к макроскопическому масшіабу сисіемьі. Выражение для дисперсии насыщенное і и можно представшь похожим образом „iJAz&L-erfifl-r&b)) (3.4.8)
Учитывая соотношение (3.4.6). можно увидеть, что выражения (3.4.7) и (3 4.8) не зависят от внешнего размерного параметра L (в нашей носіановке L — ос)
Используя последние формулы можно ючно вычислить и посіроиіь графически распределение средней насыщенное і и и дисперсию насыщенное і и вблизи фроніа. Для иросюїьі примем, чю насыщенность спереди фронa St = 0 (изначально резервуар был заполнен только вытесняемой жидкостью). В случае квадратичных функций oi носи тельных проницаемое тей (2.1.13) можно вычислить значение насыщенносш на, задней сюроне поверх-носіи фронт (см. приложение АЛ) и принять S\ равным эюму значению. с - 1 Ьі - . , л/1 + т где т = Цо/l w Для производной о г функции распределения по і оков имеем (см. приложение А.1)
При заданных функциях относительной проницаемости, значениях вязкостей вытесняющей жидкости JJLW и вытесняемой жидкости ц0 можно подсчитать F (Si) и значение константы А (см. (2 3.16)). Это дает возможность построить графически (S(r,t)) и а для требуемого значения безразмерною времени tj при заданных параметрах поля проницаемости Ко и а.
На рис 8 графически представлен профиль средней насыщенности для различных отношений вязкостей вытесняемой и закачиваемой жидкостей для одною и юю же значения безразмерного времени tj = 0,5. Как видно из рисунка, при увеличении отношения вязкостей m = \i0f\iw профиль насыщенносш оказьіваеіся более размытым, чю сооївеїсівуег большим возмущениям фронта вьпеснения. Отметим іакже. чю насыщенность на задней стороне поверхности фронта S\ в этом случае меньше, чю говорит о менее эффективном вытеснении. На рис 9 представлено сравнение дисперсии насыщенное і и для различных значений т: очевидно, что чем выше т тем больше дисперсия насьіщенносіи. На основе этих наблюдений можно сделать вывод, что на практике, в нефтяной индусірии, увеличение отношения вязкостей приведет к большим остаточным насыщенностям нефти в резервуаре и бо лее раннему появлению воды в добывающих скважинах, что окажется менее благоприятным сценарием разработки месторождения
На рис. 10 графически представлен профиль средней насыщенности для различных значений дисперсии логарифма проницаемости для одною и того же значения безразмерною времени td = 0,5. Как видно из рисунка, увеличение дисперсии логарифма проницаемосіи веде г к растяжению профиля средней насыщенности, т.е. к большим возмущениям фронт вытеснения. На иракіике. в нефіяной индустрии, вытеснение в средах с большими значениями дисперсии лоїарифма проницаемости будет менее зффекіивньїм и буде і" сопровождайся более ранним появлением воды в добывающих скважинах. На рис. 11 представлено сравнение дисперсии насыщенноаи для различных значений Ко. Чем выше KQ іем больше дисперсия насыщенности.
Дискретизация члена, учитывающего неоднородное в среды
Для тою. чюбы вычислип» член Vb в сооїветствии с (4.2.4). необходимо перейти в пространство Фурье-образов, умножить Фурье-образ функции положения фронта на модуль волнового вектора и затем вернуться в реальное пространство. Для осуществления перехода в пространство Фурье-образов используется дискретное преобразование Фурье, которое может бьиь выполнено по определенному алгоритму ускоряющему вычисления и называемому «бьісірое преобразование Фурье» (наибольшей эффекхивносш алгоритм быстрого преобразования Фурье достигает, если число ячеек, по которым производиіся преобразование, равняется 2П. где п — целое положительное число, см. [86. 76]).
Двумерное дискретное преобразование Фурье функции Sh по поперечным координаїам можно определить следующим образом Sh{qil43) = 6-2 -1)( -1)/ ,(,,-1)0-1)/ (43.з) k-l j=l где qi и qs в дискреіном пространстве соответствуют компонентам волновою вектора q в непрерывном пространстве, a j и к соотве їству ют поперечным координатам у и z в непрерывном нроорансгве. Для простоты рассмотрим случай, когда Ny и Nz четные (процедура распросіраняеіея с незначительными изменениями на случай неченюю Ny или Nz). Учитывая симметричный порядок, в коюром хранятся данные при использовании быстрого преобразования Фурье (см. [87J). умножению на модуль волновою вектора в непрерывном пространстве (см. (4.2.4)) сооївеїсгвует следующий алюріпм в дискретом просі ранстве Фурье-образов (6hq обозначает функцию положения фроша умноженную на модуль волновою вектора): где оператор цикла for действует так же, как в языках программирования (происходит увеличение индекса на единицу). Таким образом «пробегается» весь диапазон значений q и д. . производя умножение 5h{qi,qi) на модуль ВОЛИОВОІ о вектора по определенному выше правилу. Тогда окончаїельно для дискретной формы записи V& получаем У 2 93-192=1
Приме]) Фурье-образа флуктуации функции положения фроша 8h{q2, уз) представлен на рис. 24 (Ny = Nz = 1G изображена вещественная часіь функции). Результат умножения на модуль волновою вектора представлен на рис. 25. Из рисунка видно, чю амплиіудьі возросли в соответствии с правилом умножения и с учеюм симметрии записи данных при использовании бьісірого преобразования Фурье. Таким образом, член V на каждом временном шаге будет ответственен за изменение формы фронта в зависимости от параметра Л (те. в зависимосіи от отношения подвижностей на фронте).
В качестве шага по времени выбираем при -пом обеспечивается прохождение фронтом каждой ячейки сеіки.
В вьічислиіельньїх экспериментах по исследованию течений в пористых средах часто используеіся так называемое безразмерное время, которое пред-ствляет собой объем закачанной жидкости в единицах объема норового пространства резервуара
По предложенному алгоритму численною метода расчеіа формы фроша была разработана программа. Предсіавление резулыаюв работы нашей про-1 раммы приводи іся в форме сравнения с резулг/гатами. даваемыми уже суіце-ствуюіцей и широко используемой программой моделирования двухфазных потоков. Была выбрана программа 3DSL [54. 55]. основанная на линиях тока, а следовательно обладающая существенно сниженной численной дисперсией. Сравнение представлено на примере неоднородной сетки проницаемое і и, имеющей размеры Nx = Ny = Nz = 32 ячейки и сгенерированной с использованием корреляционной функции Гауссовой формы с дисперсией логарифма проницаемости KQ = 0,1 И длиной корреляции а = б ячеек. В программе 3DSL проводился полный численный эксперимент но двухфазному вытеснению для заданною отношения вязкостей т с использованием квадратичных функций относительной проницаемости (2.1.13). Для достижения достаточно высокой іочное і и вычислений в каждом вычислительном эксперименте с использованием 3DSL иоле давления обновлялось 20 раз (отметим, что точность результатов, даваемых программой 3DSL может сильно зависеть от числа обновлений поля давления). Затем, для значения безразмерного времени lci = 0,49 определялось положение фронта, которое изображалось іра-фически. После этого, для той же реализации поля проницаемости, форма фроша определялась с помощью предлагаемой программы быстрого моделирования для гого же значения td = 0,49 и отношения вязкосіей т и также изображалась графически. Затем оіношение вязкостей т изменялось, вычислительные эксперименты в обеих программах повторялись и снова проводилось сравнение полученных результатов. Результаты сравнения попарно преде ывлены на рисунках, приводимых ниже.
На рис 26 и 27 представлено сравнение для т = 0,5. когда вязкость закачиваемой жидкости в два раза больше вязкости вытесняемой жидкости. Как видно из рисунка, фронт ночі и плоский, с незначиїельньши возмущениями. Результат рабоїьі предлагаемой программы, хорошо соответствует результату, получаемому с помощью программы 3DSL.
На следующих рис. 28 и 29 предсіавленьї результаты после изменения отношения вязкостей таким образом, что вязкость закачиваемой жидкости сіала в два раза меньше вязкости вытесняемой жидкости (гп = 2). Как видно из рисунка, фронт, полученный с помощью предлагаемой программы моделирования, изменился в соответствии с изменением отношения вязкостей (амплитуда возмущений увеличилась). Это хорошо соотносится с результатом, полученным после выполнения полного численного эксперимента в программе 3DSL
И наконец, на рис. 30 и 31 приводятся результаты работы программ при отношении вязкостей т = 3 которому соответствует отношение подвижно-стей на фронте М/ = 1 при использовании квадратичных функций относительной проницаемости. Такое отношение параметров является критическим при рассмотрении устойчивого фронта. Как видно из рисунков, соответствие результатов, даваемых двумя программами, хорошее.
