Solution of the problem with a small parameter in suspension filtration in a porous medium

Texto integral

Введение

Большинство физических процессов описывается линейными дифференциальными уравнениями в частных производных первого и второго порядка [1–4]. Рассматриваемые процессы протекают в пространстве и времени, а уравнения содержат соответствующие производные. Поэтому для однозначности решения ставится задача Коши, заключающаяся в определении состояния системы по её известной начальной конфигурации [1, 5–7], а также краевая задача на границах рассматриваемой области [6, 8, 9].

Уравнения параболического типа описывают монотонное распространение физического параметра в рассматриваемой области. В частности, к уравнениям параболического типа [6] относятся уравнения массо- и теплопереноса – уравнения диффузии и теплопроводности [2, 10, 11]. В ряде задач наблюдаются разрывные распределения физических величин в области моделирования [12–14]. Такое поведение решения справедливо для уравнений гиперболического типа первого порядка, в частности для закона сохранения массы [3, 15].

Решения дифференциальных уравнений в частных производных первого и второго порядка для описания физических процессов являются корректными в случае, если сами уравнения являются эволюционными, то есть дающими истолкование этих процессов как последовательную смену некоторых состояний без наличия существенных флуктуаций [16]. В случае системы дифференциальных уравнений в частных производных их эволюционность определяется из решения задачи Штурма–Лиувилля вещественностью собственных значений матриц коэффициентов [17].

В ряде случаев при описании физических систем необходимо учитывать процессы, происходящие после эволюционной смены состояний. В таком случае задача сводится к интегро-дифференциальным уравнениям типа Вольтерра или к дифференциальным уравнениям с запаздывающим аргументом [18].

Большинство физических задач в общей постановке не имеет аналитических решений, либо поиск таких решений требует больших вычислительных мощностей и временных затрат. Преимуществом эволюционных уравнений является их пригодность к использованию численных методов решения, так как в случае малого шага расчётной численной сетки по времени описывается последовательная смена близких состояний физической системы, определяемых начальными условиями [19]. Такие малые шаги по времени можно рассматривать как малые параметры при решении эволюционных систем уравнений.

Актуальность решения задач с малым параметром обусловлена разномасштабностью возникающих процессов в механике многофазных систем, например, при малообъёмных закачках физико-химических реагентов для управления процессом заводнения в нефтяном пласте [20]. В таких задачах в качестве малого параметра выступает размер оторочки реагента [20]. Малый параметр может стоять как коэффициент перед старшей производной (такие задачи называются сингулярно возмущёнными) или в граничных условиях.

С помощью условий, накладываемых на значения малого параметра, для рассматриваемого вида уравнений необходимо контролировать выполнение условий теоремы Коши о существовании и единственности решения [21].

Для решения задач с малым параметром необходимо использование специальных методов теории дифференциальных уравнений, аппарата математического анализа и математической физики [22], в том числе учитывающих возможность эффекта пограничного слоя [23]. При наличии такого эффекта определённые значения малого параметра могут приводить к существенным изменениям решения вблизи границы рассматриваемой области [23], что и наблюдается в процессе применения физико-химических методов увеличения нефтеотдачи пластов [20], в частности суспензий. Кроме того, решение таких уравнений непрерывно зависит от значения малого параметра [24, 25].

Распространённым подходом для моделирования задач подземной гидромеханики является использование гидродинамических симуляторов [26, 27]. В этих симуляторах хорошо себя зарекомендовали методы механики сплошных сред. Однако в настоящее время не исследована корректность применения этих симуляторов для решения задач с малым параметром.

Другим распространённым подходом к моделированию процессов фильтрации суспензий в пласте является использование модели глубокого проникновения частиц реагента в пористую среду [28–30]. В этой модели снижение пористости и проницаемости пласта из-за оседания дисперсных частиц описывается введением двух параметров – коэффициентов фильтрации и повреждения породы. На практике эта модель используется для описания процессов в ходе применения технологии выравнивания профиля приёмистости (выравнивания скоростей фильтрации и расходов по пропласткам) нагнетательной скважины. Целесообразным видится её применение также для описания схожего процесса кольматации пористой среды при заводнении. Глубина проникновения дисперсных частиц при этом не превышает нескольких метров [20], следовательно, эта глубина и является малым параметром, что отмечается в данной работе впервые. Целью настоящей работы является решение задачи фильтрации суспензии в пористой среде в локальной области вблизи добывающей скважины. Впервые в рамках единой модели рассматриваются процессы выравнивания профиля приёмистости нагнетательной скважины и кольматации при заводнении, указываются малые параметры в каждом случае.

Математическая модель глубокого проникновения дисперсных частиц в пористую среду

Представлена модель однородной пористой среды, насыщенной раствором полимера (рис. 1).

 

Рис. 1. Пористая среда. Обозначения: 1 – порода; 2 – раствор полимера; 3 – частицы в потоке; 4 – частицы, удержанные на скелете пористой среды

Fig. 1. Porous medium. Description: 1 – rock; 2 – polymer solution; 3 – particles in flow; 4 – particles retained on the rock matrix

 

Концентрация дисперсных частиц в растворе полимера обозначается массовой концентрацией c. В то же время часть частиц σ, относящаяся к объему пористой среды, задерживается в сужениях поровых каналов (на скелете пористой среды).

Предполагается, что скорость раствора полимера и взвешенных в нем частиц является одинаковой и равной u (частицы полностью увлекаются несущей фазой), а также если пренебречь диффузионной частью в выражении для плотности потока, что справедливо для достаточно крупных частиц [20], закон сохранения массы частиц записывается в виде:

$\frac{\partial \left(\left(m_0-\sigma \right)c+\sigma \right)}{\partial t}+\overrightarrow{\nabla }\left(\overrightarrow{u}c\right)=0$, (1)

где t – время; m₀ – начальная пористость.

Первое слагаемое в (1) учитывает тот факт, что при осаждении частиц значение пористости уменьшается [20], и определяет, сколько в единице объема пористой среды содержится частиц.

Скорость осаждения частиц зависит от концентрации частиц на единицу объема пористой среды, массовой концентрации частиц в полимере и скорости фильтрации [20]. Обычно полагается, что полимер присутствует во всем потоке с постоянной начальной концентрацией. В рассматриваемом процессе полимерная составляющая играет роль «цемента», который закрепляет захваченные частицы на скелете породы. С другой стороны, полимер удерживает свободные частицы во взвешенном состоянии. В результате в скорости захвата частиц не будет члена, отвечающего за срывание частиц и затягивание их обратно в поток [20].

Для описания динамики концентрации удержанных частиц можно использовать классическое уравнение глубокого проникновения частиц [20, 28–30]:

$\frac{\partial \sigma }{\partial t}=\lambda uc$, (2)

где λ – коэффициент фильтрации.

Коэффициент фильтрации в (2) определяют экспериментально [30] и нередко считают постоянным. Однако в общем виде λ считают функцией от концентрации осажденных частиц в единице объема пористой среды [28, 29]:

$\lambda ={\lambda }_0\left(1-\frac{\sigma }{m_0}\right)$, (3)

где λ₀ – коэффициент фильтрации до закачки суспензии полимера.

Для описания течения раствора полимера в пористой среде можно использовать модифицированное уравнение Дарси, учитывающее снижение проницаемости за счет осаждения частиц [20]:

$\overrightarrow{u}=-\frac{k}{\mu \left(1+\beta \sigma \right)}\text{grad}p$, (4)

где k – абсолютная проницаемость пористой среды; µ – коэффициент динамической вязкости; p – поровое давление; β – коэффициент повреждения породы. Коэффициент β определяется геометрией поровой среды и показывает степень загрязнения пористой среды [20].

Таким образом, с учётом (1), (3), (4) в общем виде система уравнений глубокого проникновения дисперсных частиц в пористую среду имеет вид:

$\frac{\partial \left(\left(m_{0i}-{\sigma }_i\right)c_i+{\sigma }_i\right)}{\partial t} +\overrightarrow{\nabla }\left({\overrightarrow{u}}_i{c}_i\right)=0$, (5)

$\frac{\partial {\sigma }_i}{\partial t}={\lambda }_0\left(1-\frac{{\sigma }_i}{m_{0i}}\right)u_i{c}_i$, (6)

${\overrightarrow{u}}_i=-\frac{k_i}{\mu \left(1+\beta {\sigma }_i\right)}\text{grad}p$, (7)

где i – номер пропластка.

Система уравнений (5)–(7) далее будет называться полной. Эта система уравнений может быть решена численно, например, при помощи явной схемы, так как аналитически такую задачу решить не получится ввиду ее нелинейности. Однако нужно понимать, что на промысле точность и время – это главные параметры, позволяющие правильно и быстро предсказывать дальнейшую работу участка месторождения. К сожалению, для повышения точности расчётов приходится выбирать достаточно малые шаги по координатам численной сетки в случае использования численных решений или гидродинамических симуляторов, что приводит к росту времени расчётов.

Если предположить, что концентрация осажденных частиц на единицу объема пористой среды много меньше начальной пористости (σ<<m₀), можно получить упрощенную систему уравнений глубокого проникновения частиц [20]:

$\frac{\partial \left(m_{0i}{c}_i+{\sigma }_i\right)}{\partial t}+\overrightarrow{\nabla }\left({\overrightarrow{u}}_i{c}_i\right)=0$, (8)

$\frac{\partial {\sigma }_i}{\partial t}={\lambda }_0{u}_i{c}_i$, (9)

${\overrightarrow{u}}_i=-\frac{k_i}{\mu \left(1+\beta {\sigma }_i\right)}\text{grad}p$. (10)

Упрощенная система уравнений м быть решена аналитически [31]. Аналитическое решение, в отличие от численного, будет давать точный и быстрый результат.

Таким образом, получены две системы уравнений: полная – учитывает влияние осаждаемых частиц на значение пористости коэффициент фильтрации и может быть решена численно, упрощенная – пренебрегает размером осадка по отношению к начальной пористости может быть решена аналитически.

Сопоставление решений полной и упрощенной систем уравнений

Для проверки предположения о малости концентрации осажденных частиц в единице объема пористой среды по сравнению с начальной пористостью необходимо решить полную и упрощенную системы уравнений – (5)–(7) и (8)–(10), соответственно.

Пусть рассматривается задача о полимер-дисперсной обработке слоисто-неоднородного пласта вблизи трещины. Эта задача подробно описана и аналитически решена при помощи метода характеристик с использованием упрощенной системы уравнений глубокого проникновения частиц для линейного (одномерного) течения в работе [20].

Поэтому разумно использовать результаты авторов [20], а также их параметры модели для того, чтобы сопоставить решение полной системы уравнений с решением упрощенной. Как выше было сказано, полная система может быть решена численно, например, при помощи явной конечно-разностной схемы.

Предварительно нужно обезразмерить систему уравнений. Для этого вводятся следующие безразмерные переменные [20]:

$X=\frac{x}{L}$, (11)

$\Lambda ={\lambda }_{0}L$, (12)

$T=\frac{Qt}{2LHl}$, (13)

$V_i=\frac{k_{i}H}{{\Sigma }_j{k}_j{h}_j}$, (14)

где L – характерный линейный размер задачи; Q – расход суспензии в нагнетательной скважине; l – длина трещины; H – эффективная толщина пласта; h – толщина пропластка; x – горизонтальная координата; V_i – безразмерная проводимость; X, Λ и T – безразмерные координата, коэффициент фильтрации и время соответственно; j – номер пропластка.

Перед обезразмериванием вводится поток в i-ом пропластке [20]:

$Q_i=\frac{k_i{h}_i}{{\displaystyle {\sum }_j{k}_j{h}_j}}Q$. (15)

Обезразмеренная с учётом (11)–(15) система уравнений будет иметь вид

$$\begin{gathered} \frac{\partial c_i}{\partial T}+\frac{V_i}{\left(m_{0i}-{\sigma }_i\right)}\frac{\partial c_i}{\partial X}=\frac{\Lambda V_i{c}_i}{m_{0i}}\left(c_i-1\right), \\ \frac{\partial {\sigma }_i}{\partial t}=\Lambda \left(1-\frac{{\sigma }_i}{m_{0i}}\right)V_i{c}_i. \end{gathered}$$

Начальные и граничные условия для системы будут иметь вид:

$c(X=0, T\ge 0)=c_0, c(X>0, T=0)=0, \sigma (X\ge 0, T=0)=0$,

где c₀ – начальная концентрация частиц в потоке.

Решение будет производиться с использованием явной конечно-разностной схемы по программе, написанной на языке программирования Python.

Расчет будет идти до тех пор, пока T<T₀, где T₀ – безразмерный объем закачки, определяющийся по формуле

$T_0=\frac{V_{\text{зак}}}{2LHl}$,

где V_зак – объем закачки.

Результатом применения технологии выравнивания профиля приёмистои является снижение проницаемости в призабойной зоне, что позволяет перераспределить потоки гнетаемой жидкости между пропластками с низкой и высокой проницаемостью, то есть провести выравнивание скоростей фильтрации в пропластках, а следовательно, и фронта вытеснения нефти водой.

Проницаемость для i-го пропластка после обработки определяется формулой:

$k_{1i}=\frac{k_i}{{\displaystyle {\int }_0^1\left(1+\beta {\sigma }_i\right)dX}}$. (16)

Знаменатель в формуле (16) называется дополнительным фильтрационным сопротивлением g_i, которое образовалось в результате закрепившихся частиц на скелете породы. Так как из-за неоднородности пласта в каждый пропласток входит различный объем реагента, дополнительное фильтрационное сопротивление для каждого пропластка будет разным.

Для того чтобы сопоставить значения дополнительного фильтрационного сопротивления, определяемого при помощи аналитического решения и численного, интеграл в формуле (16) можно решить, например, при помощи метода Монте-Карло.

Пусть имеется модельный п с параметрами, приведенными в таблице. Всего имелось пять пропластков. Значения эмпирических параметров были следующими: λ₀=10 м^–1, β=400. Начальная концентрация суспензии взята c0,02, а для вязкости принималось значение μ=0,8 Па∙с. Длина трещины l=79 м, а размер призабойной зоны скважины L=40 м. Объем закачки брался равным 1302 м³ [20]. Шаг по координате ΔX=10^–3 и по времени ΔT=10^–5.

 

Таблица. Параметры пропластков

Table. Interlayers pameters

Номер пропластка/Interlayer no.	ki, мД/mD	hi, м/m	m0i
1	1213	1,3	0,236
2	371	1,4	0,209
3	52	4,7	0,189
4	106	1,6	0,201
5	135	1,3	0,196

                                              

Рассматривалась обработка вертикальной скважины с трещиной гидроразрыва. Распределения частиц в потоке и доли осажденных частиц на скелете породы в наиболее высокопроницаемом пропластке с проницаемостью 1213 мД приведены на рис. 2, 3 соответственно. Для сравнения также приведены распределения концентраций частиц вблизи трещины в высокопроницаемом пропластке, полученные аналитически с использованием упрощенной системы уравнений [20].

 

Рис. 2. Распределение концентрации взвешенных частиц в потоке вблизи трещины. Решения полной (синяя кривая) и упрощенной (красная кривая) систем уравнений

Fig. 2. Distribution of volumetric particle content moving with the flow near a fracture. Solutions of complete (blue curve) and simplified (red curve) systems of equations

 

Рис. 3. Распределение концентрации осажденных частиц вблизи трещины. Решения полной (синяя кривая) и упрощенной (красная кривая) систем уравнений

Fig. 3. Distribution of the volumetric particle content retained in the unit volume of porous media near a fracture. Solutions of complete (blue curve) and simplified (red curve) systems of equations

 

Как видно из рисунков, решения, полученные для полной и упрощенной систем уравнений, соответствуют друг другу с относительной ошибкой менее 4 %. Действительно, концентрация осажденных частиц в единице объема пористой среды мала по сравнению с начальной пористостью. Этот результат позволяет использовать данное предположение для упрощения системы уравнений глубокого проникновения частиц, которую в дальнейшем можно решить аналитически.

Также были сопоставлены значения дополнительного фильтрационного сопротивления, рассчитанного аналитически при помощи упрощенной системы уравнений и рассчитанного численного при помощи полной.

Для высокопроницаемого пропластка дополнительное фильтрационное сопротивление, полученное при помощи решения упрощенной и полной систем уравнений, равно 1,744 и 1,807 соответственно. Разница между значения примерно равна 3,5 %. В некоторых ситуациях такое отклонение может быть значительным, поэтому был проведен расчет для шага по координате ΔX=10^–4 и по времени ΔT=10^–6 для того, чтобы численное решение стремилось к аналитическому. В таком случае значения дополнительного фильтрационного параметра равны 1,744 и 1,740 соответственно для упрощенной и полной систем уравнений. Теперь разница ­между значениями составляет 0,2 %. Данный результат также подтверждает допущение о малости концентрации осажденных частиц по сравнению с начальной пористостью и показывает корректность полученного решения полной системы уравнений при помощи явной конечно-разностной схемы (при уменьшении шагов численной сетки решение должно стремиться к аналитическому – условие сходимости).

Математическая модель кольматации пористой среды

В другом подходе учитывается изменение пористости вследствие осаждения частиц в поровых каналах. В таком случае скорость осаждения частиц и снижения пористости характеризуется коэффициентом кольматации, который также является малым параметром. Система уравнений при этом остаётся подобной уравнениям глубокого проникновения дисперсных чтиц в пористую среду, но меняется закон осаждения частиц и соответствующее уменьшение пористости [32].

Фильтрация суспензий с небольшой массовой долей дисперсных частиц описывается системой уравнений, в которую входит баланс массы частиц в суспензии и несущей жидкости, закон Дарси, функциональная зависимость проницаемости от пористости, например, рмула Козени–Кармана и уравнение, описывающее динамику оседания частиц на скелете пористой среды [32]:

$\begin{array}{c}\frac{\partial }{\partial t}\left(mc\right)+div\left(c\overrightarrow{u}\right)=\frac{\partial m}{\partial t} ,\\ \frac{\partial }{\partial t}\left(\left(1-c\right)m\right)+div\left(\left(1-c\right)\overrightarrow{u}\right)=\mathrm{0,} \\ \overrightarrow{u}=-\frac{k\left(m\right)}{\mu }gradp , k\left(m\right)=k_0{\left(\frac{m}{m_0}\right)}^r,\\ \frac{\partial m}{\partial t}=-\gamma c\left(m-m_{\text{ст}}\right),\end{array}$

где с – объемная концентрация частиц в суспензии; μ – вязкость суспензии; m – пористость; m_ст – стационарное значение пористости, меньше которой пористость не снижается; γ – коэффициент кольматации; k₀ – начальное значение проницаемости пласта; r – показатель степени. Предполагается, что отложение частиц на скелете не приводит к задерживанию жидкости в порах (как следствие, отсутствуют капиллярные силы), а скорость частиц совпадает со скоростью жидкости, поскольку концентрация частиц с мала.

Вводится комплекс безразмерных параметров и переменных: (11),

$T=\frac{t}{t_0}$,

$P=\frac{p-p_l}{p_w-p_l}$,

$\Gamma =\frac{\gamma }{{\gamma }_0}$,

$U=\frac{u}{u_0}$,

где P – безразмерное давление; Γ – безразмерный коэффициент кольматации; U − безразмерная скорость; p_w – забойное давление; p_l – давление на контуре питания; u₀ – характерная скорость; t₀ – характерное время; γ₀ – характерный коэффициент кольматации.

После их подстановки в исходную систему получается удобным ввести характерные значения параметров как:

$u_0=\frac{k_0\left(p_w-p_l\right)}{\mu L}$,

$t_0=\frac{\mu L^2}{k_{0(}{p}_w-p_l)}$,

${\gamma }_0=\frac{1}{t_0}$.

Тогда система уравнений в безразмерном виде будет выглядеть так:

$\frac{\partial \left(mc\right)}{\partial T}+\frac{\partial \left(Uc\right)}{\partial X}=\frac{\partial m}{\partial T}, \frac{\partial \left(\left(1-c\right)m\right)}{\partial T}+\frac{\partial \left(\left(1-c\right)U\right)}{\partial X}=0$, (17)

$\frac{\partial m}{\partial T}=-\Gamma c\left(m-m_{\text{ст}}\right)$.

Решение будет производиться с использованием явной конечно-разностной схемы.

Чтобы определить распределения объемной концентрации частиц в суспензии, необходимо учесть условия на границе пласта и начальный момент времени. Начальное условие принимается следующим:

$c\left(X\mathrm{,0}\right)=0$.

Граничное условие записывается следующим образом:

$c\left(\mathrm{0,}T\right)=c_0$

Нахождение распределения пористости также требует задания начального условия:

$m\left(X\mathrm{,0}\right)=m_0$.

В итоге имеется замкнутая система уравнений, позволяющая определить распределения пористости и концентрации дисперсных частиц в пласте.

Результаты исследования

В результате решения задачи были получены графические зависимости пористости призабойной зоны, объемной концентрации частиц в суспензии, проницаемости и давления от координаты в различные моменты времени: t=30 сут, t=60 сут, t=90 сут, остальные параметры соответствуют моделям задач подземной гидромеханики: L=100 м, µ=1 мПа·с, m_ст=0,05, γ=0,01 1/с, k₀=100 мД, r=3.

Из рис. 4 можно увидеть, что за фронтом вытеснения концентрация частиц становится равной концентрации частиц в закачиваемой смеси, а пористость уменьшается до минимального значения, которое называется стационарной пористостью. Это происходит из-за оседания дисперсных частиц и соответствующего снижения фильтрационно-емкостных свойств пласта.

 

Рис. 4. Распределение пористости и объемной концентрации частиц примеси по координате в различные моменты времени

Fig. 4. Distribution of porosity and volumetric particles concentration along the coordinate at different moments of time

 

Помимо этого, можно заметить, что фронт вытеснения продвигается дальше со временем, причём его скорость движения остается постоянной.

Также был произведен расчет распределения пористости по координате при различных значениях коэффициента кольматации $\text{Γ}: \mathrm{0,2}*{10}^3, 2*{10}^3, 20*{10}^3, 200*{10}^3$ (рис. 5), остальные параметры соответствуют модельному случаю.

 

Рис. 5. Распределение пористости по координате при различных значениях коэффициента кольматации в момент времени 30 сут

Fig. 5. Porosity distribution along the coordinate at different values ​​of the colmatation coefficient at the time of 30 days

 

При уменьшении коэффициента кольматации на фронте вытеснения не возникает разрывное решение, а образуется стабилизированная зона. В этой зоне концентрация постепенно уменьшается до нуля, а пористость увеличивается. Ширина стабилизированной зоны определяется между точками с наименьшей и наибольшей пористостью. Путем анализа результатов было выяснено, что причиной образования стабилизированной зоны является конечная скорость кольматации, которая характеризуется коэффициентом кольматации. Согласно рис. 5, чем меньше этот коэффициент, тем шире становится стабилизированная зона, что обусловлено уменьшением пористости согласно (17). Это приводит к постепенному снижению пористости по координате. Таким образом, коэффициент кольматации является малым параметром, который значительно влияет на процесс вытеснения и определяет его характер. Аналогом такого малого параметра являются капиллярные силы в задаче Рапопорта-Лиса [33].

Выводы

1. Показано, что процессы выравнивания профиля приёмистости и кольматации пористой среды описываются в рамках единой системы уравнений глубокого проникновения дисперсных частиц в пористую среду.
2. Установлено, что концентрация осаждённых частиц является малым параметром, позволяющим свести полную систему уравнений глубокого проникновения дисперсных частиц в пористую среду к упрощённому виду, при котором решение может быть получено аналитически с помощью метода характеристик.
3. Проведено сопоставление решений полной и упрощённой систем уравнений глубокого проникновения дисперсных частиц в пористую среду, свидетельствующее об их соответствии с погрешностью менее 4 %.
4. Получены распределения концентрации дисперсных частиц и пористости в пласте при кольматации, показывающие, что с течением времени фронт кольматации распространяется с постоянной скоростью.
5. Показано, что коэффициент кольматации является малым параметром, существенно определяющим характер вытеснения нефти водой, причём уменьшение коэффициента кольматации приводит к уменьшению скорости вытеснения и увеличению размеров стабилизированной зоны вблизи фронта вытеснения.

Sobre autores

Ekaterina Arkhipova

University of Tyumen

Autor responsável pela correspondência
Email: stud0000222738@study.utmn.ru
ORCID ID: 0009-0003-0022-5605

Student

Rússia, Tyumen

Alexander Gilmanov

University of Tyumen

Email: a.y.gilmanov@utmn.ru
ORCID ID: 0000-0002-7115-1629

Cand. Sc., Associate Professor

Rússia, Tyumen

Alexander Shevelev

University of Tyumen

Email: a.p.shevelev@utmn.ru
ORCID ID: 0000-0003-0017-4871

Cand. Sc., Professor

Rússia, Tyumen

Bibliografia

Arquivos suplementares