Моделирование фильтрации и рассеянного разрушения в трещиноватых зонах с аномально высоким пластовым давлением

Полный текст

Введение и обзор

Доля углеводородов, извлекаемых из нетрадиционных месторождений, таких как горючие сланцы, увеличивается из года в год. Нефтематеринские породы этих месторождений характеризуются сверхнизкой проницаемостью и аномально высоким пластовым давлением (АВПД). Залогом успешной добычи углеводородов из таких пород является наличие в окрестности добывающей скважины системы гидродинамически связанных трещин техногенного или природного происхождения. В настоящей работе рассматривается эффект от дополнительного разрушения пород пласта, связанного с АВПД.

Аномально высокое пластовое давление – давление в порах горной породы, превышающее гидростатическое давление. АВПД достаточно распространенное явление в пластах, сложенных различными породами. В нефтематеринских или газоматеринских породах АВПД формируется в основном за счет флюидогенеза (процессов превращения керогена в жидкие и газообразные углеводороды). АВПД может вызвать катастрофические явления, такие как обрушение стенок скважины при бурении с недостаточно плотным буровым раствором [Ma, Holditch, 2015; Li et al., 2012]. С другой стороны, АВПД может поддерживать приток жидкости в условиях низкой проницаемости, когда общей депрессии недостаточно [Sonnenberg, Pramudito, 2009]. Продуктивность скважин на сланцевых месторождениях, пробуренных в сходных условиях, может проявлять выраженную изменчивость от скважины к скважине. Скважины, попавшие в трещиноватую зону с АВПД, имеют высокие начальные дебиты, в то время как остальные остаются практически сухими [Алексеев, 2013; Baihly et al., 2010; Панарин, Фомин, 2016]. Кроме того, в случае остановки скважины в рамках борьбы с неконтролируемыми проявлениями АВПД при повторном запуске не удается добиться прежних высоких дебитов [Алексеев, 2013]. Все это свидетельствует о том, что в районе скважин идут необратимые процессы разрушения, связанные с наличием АВПД. Таким образом, АВПД можно рассматривать как потенциально полезный дополнительный источник пластовой энергии, освоение которого является актуальной задачей.

В частности, большой интерес представляет влияние рассеянного трещинообразования на дебит при добыче углеводородов из низкопроницаемых пластов с АВПД. Эффект растрескивания матрицы под действием АВПД был описан в литературе достаточно давно. В работе [Secor, 1965] явлению саморастрескивания под действием АВПД в непроницаемой среде было дано название автофлюидоразрыв (natural hydraulic fracturing, NHF). Позже было сделано обобщение на случай пороупругой среды [Engelder, Lacazette, 1990]. Утверждается, что это явление важно для процессов миграции жидкости на значительных (геологических) масштабах времени. Авторазрыв возникает, когда постепенно увеличивающееся давление жидкости в порах достигает некоторого критического предела [Ougier-Simonin et al., 2016]. Очевидно, что аналогичный эффект может быть достигнут во время искусственного воздействия на пласт с АВПД, когда начальное напряженное состояние в нем тем или иным образом нарушается (бурение, гидравлическое разрыв и т.д.). В настоящей работе рассматривается случай, когда АВПД уже сформировалось, но критическое давление для развития авторазрыва еще не достигнуто.

Теория повреждаемости (CDM) является естественным подходом к описанию зарождения, роста и слияния микродефектов, включая трещины [Lemaitre, 1996; Murakami, 2012]. При условии достаточно большого количества включений и микродефектов в каждом элементарном объеме среды описание процессов возможно в рамках предположения о сплошной среде. В теории повреждаемости в набор термодинамических параметров сплошной среды добавляется скалярный или тензорный параметр поврежденности. Система определяющих отношений включает в себя специальное кинетическое уравнения для этого параметра. Параметр поврежденности влияет на поведение сплошной среды. Развитие микродефектов приводит к деградации модулей упругости или появлению остаточных деформаций. В работе [Кондауров, Фортов, 2002] отмечается особая роль энергетического баланса в эволюции микродефектов. Впоследствии подход Кондаурова был обобщен на случай насыщенных пористых среда с хрупким скелетом [Извеков, Кондауров, 2009] и на случай насыщенной среды с двойной пористостью с хрупким скелетом [Извеков и др., 2020].

Основной моделью для континуального описания потока жидкости в трещиновато-пористой среде является модель двойной пористости, впервые предложенная в работах [Баренблатт и др., 1960; Warren, Root, 1963]. Семейство сплошных моделей среды с двойной пористостью основано на представлении среды как суперпозиции трех континуумов – твердого скелета и двух флюидов (жидкость в трещинах и жидкость в матрице). Между матрицей и трещинами возможен массообмен. Существует множество модификаций модели. Учитывается нелинейный закон массообмена [Zimmerman, 1993], деформируемость скелета [Bai, 1995; Berryman, 2002], иерархическое строение пустотного пространства скелета (множественная пористость) [Mehrabian, Abousleiman, 2015], и т.д.

В работе [Извеков и др., 2020] была представлена континуальная модель среды с двойной пористостью с повреждаемым скелетом и АВПД. Определяющие соотношения и уравнение эволюции для параметра поврежденности согласуются с фундаментальными законами сохранения и вторым законом термодинамики. Поскольку характерные напряжения в пласте много меньше модулей упругости, принимается приближение малых деформаций. Это предположение позволяет сформулировать линеаризованную систему определяющих соотношений. Однако закон массообмена остается существенно нелинейным по отношению к параметру поврежденности.

Настоящая работа является развитием [Извеков и др., 2020]. Сформулирована и численно решена начально-краевая задача в одномерной постановке, моделирующая связанные процессы течения, разрушения и изменения напряженно-деформированного состояния в пласте, изначально находившемся в равновесии при АВПД. Постановка задачи соответствует процессам в окрестности плоской расклиненной трещины гидроразрыва на стадии добычи. Исследованы различные режимы распространения волн разрушения в пласте. Учтена конечная протяженность зоны проводящих трещин. Отметим, что в общем случае сланцы представляют собой анизотропную среду. Более того, модули упругости, прочностные свойства и проницаемость демонстрируют анизотропные свойства [Ougier-Simonin et al., 2016; Dokhani et al., 2016]. Однако учет анизотропии не влияет на основные особенности процесса распространения волн разрушения. Поэтому в настоящей работе для простоты среда считается изотропной.

Система уравнений

Подробно построение модели рассматриваемой среды и оценка входящих в нее констант описана в работе [Извеков и др., 2020]. Приведем здесь основные результаты.

Рассматривается насыщенная одной жидкостью среда с двойной пористостью, матрица (блоки) которой способна накапливать поврежденность. Среда моделируется совокупностью трех взаимопроникающих континуумов. Две жидкие сплошные среды представляют собой результат усреднения жидкости, насыщающей блоки и трещины, третий представляет собой скелет – сплошное деформируемое повреждаемое твердое тело. Жидкие среды могут обмениваться массой. Скелет в массообмене не участвует. Предполагается, что поврежденность, накапливаемая в блоках матрицы, в первую очередь влияет на закон массообмена. Рассматриваемая модель линеаризуется, однако обменный член остается принципиально нелинейным по отношению к параметру поврежденности.

Система уравнений рассматриваемой среды состоит из уравнений неразрывности для континуумов (1), (2), уравнения совместности деформаций и скоростей для скелета (3), законов движения флюидов (здесь в форме закона Дарси) (4) и уравнения равновесия среды как целого (5):

$\frac{\partial r_1}{\partial t}+\nabla \cdot \left(r_1{v}_1\right)={\dot{m}}_{12},$     (1)

$\frac{\partial r_2}{\partial t}+\nabla \cdot \left(r_2{v}_2\right)={\dot{m}}_{21},$     (2)

$\frac{\partial e}{\partial t}=\frac{1}{2}\left(\nabla \otimes v_S^{}+\nabla \otimes v_S^T\right)$     (3)

$W_{\alpha }=-\frac{k_{\alpha }}{{\mu }_f}\left(\nabla p_{\alpha }-{\rho }_{\alpha }g\right),$     (4)

$\nabla \cdot T=-rg,$     (5)

где: индекс α означает номер флюида (1 – в матрице, 2 – в системе трещин); ${\phi }_{\alpha }$ – объемная доля флюидов; $r_{\alpha }={\phi }_{\alpha }{\rho }_f$ – плотность флюидов; ${\rho }_f$ – плотность жидкости; $v_{\alpha },v_S$ – скорости флюидов и скелета; ${\dot{m}}_{12}=-{\dot{m}}_{21}\equiv \dot{m}$ – скорость массообмена между матрицей и трещинами в единице объема, скорость фильтрации $W_{\alpha }={\phi }_{\alpha }{w}_{\alpha }={\phi }_{\alpha }\left(v_{\alpha }-v_s\right);$ $p_{\alpha }$ – давление жидкости в подсистемах; $k_{\alpha }$ – проницаемость; ${\mu }_f$ – вязкость жидкости; g – ускорение свободного падения; $T=T_s-\left(p_1{\phi }_1+p_2{\phi }_2\right)I$ – тензор Коши полных напряжений в пористой среде; $T_s$ – тензор парциальных напряжений в скелете; $T_{\alpha }=-p_{\alpha }{\phi }_{\alpha }I$ – тензор парциальных напряжений во флюидах; e – тензор малой деформации скелета; $r=\sum _A{r}_A,$ $I={\delta }_{ij}{e}_i\otimes e_j$  – единичный тензор второго ранга; ${\delta }_{ij}$ – символ Кронекера; ei – векторы ортонормированного базиса, знаком “⊗” обозначено тензорное умножение.

Приведенное неравенство Клаузиуса–Дюгема (второе начало термодинамики) для рассматриваемой среды в изотермическом приближении имеет вид [Извеков и др., 2020]:

$\begin{array}{c}-\sum _A{r}_A\frac{d_A{\psi }_A}{dt}+T_s:\frac{d_{s}e}{dt}+\\ +\sum _{\alpha }\left(p_{\alpha }\frac{d_s{\phi }_{\alpha }}{dt}+\frac{p_{\alpha }{\phi }_{\alpha }}{{\rho }_{\alpha }}\frac{d_{\alpha }{\rho }_{\alpha }}{dt}\right)+{\delta }_w+{\delta }_m\ge 0.\end{array}$     (6)

Здесь введены обозначения: ${\psi }_A=u_A-{\eta }_A\theta$ – удельная свободная энергия континуума A; $u_A,$ ${\eta }_A$ – удельные внутренняя энергия и энтропия; температура $\theta$ в подсистемах считается одинаковой; ${\delta }_w=\sum _{\alpha }\frac{{\mu }_{\alpha }}{{\phi }_{\alpha }{k}_{\alpha }}{\left|W_{\alpha }\right|}^2\ge 0$ – диссипация, связанная с фильтрацией флюидов, подчиняющейся закону Дарси; ${\delta }_m=-\sum _{\alpha \ne \beta }{\chi }_{\alpha }{\dot{m}}_{\alpha \beta }$ – диссипация, связанная с обменом массой между подсистемами; ${\chi }_{\alpha }={\psi }_{\alpha }+p_{\alpha }/{\rho }_{\alpha }$ – удельный потенциал Гиббса флюидов. Неотрицательность δ_m обеспечивается дополнительным условием, связывающим обмен массой с разностью удельных потенциалов Гиббса, в изотермическом случае с разностью давлений.

Система (1)–(5) замыкается системой определяющих соотношений. Согласно работам [Трусделл, 1975; Truesdell, Noll, 2004] определяющие соотношения должны удовлетворять ряду принципов, например, требуется инвариантность относительно смены системы отсчета, локальность, термодинамическая согласованность. В качестве определяющих соотношений рассматриваемой среды примем набор функций, связывающих реакцию $\Upsilon$ среды с ее текущим состоянием [Извеков и др., 2020]:

$\Upsilon =\Upsilon \left(e,{\rho }_{\alpha },w_{\alpha },\omega \right),$     (7)

где $\Upsilon =\left\{{\psi }_{\alpha },T_S,p_{\alpha },{\phi }_{\alpha },{\dot{m}}_{12}\right\}.$ Соотношения (7) должны быть дополнены кинетическим уравнением для параметра поврежденности:

$\frac{d_S\omega }{dt}=\Omega \left(e,{\rho }_{\alpha },w_{\alpha },\omega \right),$     (8)

где $\Omega$ – функция текущего состояния рассматриваемой среды.

Определяющие соотношения, необходимые и достаточные для выполнения (6) в произвольном процессе имеют следующий вид [Извеков и др., 2020]:

${\psi }_{\alpha }={\psi }_{\alpha }\left({\rho }_{\alpha }\right),$ $p_{\alpha }={\rho }_{\alpha }^2\frac{\partial {\psi }_{\alpha }}{\partial {\rho }_{\alpha }},$     (9)

$T=r_s\frac{\partial \Psi }{\partial e},$                             (10)

${\phi }_{\alpha }=-r_s\frac{\partial \Psi }{\partial p_{\alpha }}.$                      (11)

С учетом соотношений (9)–(11) неравенство (6) принимает вид:

${\delta }_w+{\delta }_m+{\delta }_{\omega }\ge \mathrm{0,}$               (12)

где

${\delta }_{\omega }=-\left(\frac{\partial \Psi }{\partial \omega }\right)\frac{d_s\omega }{dt}=-\left(\frac{\partial \Psi }{\partial \omega }\right)\Omega$    (13)

– диссипация, связанная с растрескиванием матрицы.

В выражении (10) использован потенциал скелета [Кондауров, 2007]

$\Psi \left(e_s,p_{\alpha },\omega \right)=F-\sum _{\alpha }\frac{p_{\alpha }}{r_s}{\phi }_{\alpha },$    (14)

где $F\left(e_s,{\phi }_{\alpha },\omega \right)={\psi }_s\left(e_s,{\phi }_{\alpha }\left(e_s,{\rho }_{\alpha },\omega \right),\omega \right).$      (15)

Процедура доказательства соотношений (9)–(12) аналогична в работах [Кондауров, 2007; Truesdell, Noll, 2004]. В общих чертах она заключается в следующем. При подстановке соотношений (7), (8) в (6) получается неравенство вида $\sum _i{f}_i\left(a_1,a_2,\mathrm{....}\right){\dot{a}}_i+\delta \ge \mathrm{0,}$ где $f_i$ – функция параметров состояния $a_1,a_2,\mathrm{...,}$ которая умножается на ${\dot{a}}_i$ – скорость изменения i-го параметра состояния. Согласно принципу термодинамической согласованности $f_i$ и ${\dot{a}}_i$ независимы, при этом ${\dot{a}}_i$ – произвольны. Таким образом, неравенство может быть выполнено только в случае $f_i\left(a_1,a_2,\mathrm{...}\right)=\mathrm{0,}$ $\delta \ge 0.$

Определяющие соотношения (9), (10) и (11) линеаризуются в окрестности особого состояния $\left\{\mathrm{0,}{\rho }_{\alpha }^0,\mathrm{0,}\omega \right\},$ в котором скелет не деформирован, жидкости покоятся, давления флюидов $p_1^0=p_2^0=0.$ Давления во флюидах и напряжения в скелете в произвольных состояниях будем считать малыми по сравнению с характерными упругими модулями скелета $p_{\alpha }/K_S<<\mathrm{1,}$ $||T_S||/K_S<<\mathrm{1,}$  где $K_S$ – модуль объемного сжатия скелета. Далее рассматривается частный случай, в котором начальные полные напряжения в состоянии $\left\{\mathrm{0,}{\rho }_{\alpha }^0,\mathrm{0,}0\right\}$ равны нулю.

Эксперименты с пиролизом керогена показывают, что повышение давления при разложении керогена и извлечении продуктов пиролиза приводит к увеличению количества и размеров микротрещин. При отсутствии внешних сжимающих напряжений наблюдается значительное увеличение пористости и объемной деформации [Saif et al., 2017; 2019]. Принимая во внимание этот факт, мы предполагаем, что основным макроскопическим эффектом рассеянного микроразрушения блоков матрицы является остаточная деформация скелета. Это допущение формально сближает данный вариант модели поврежденности с теорией пластичности.

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

$T=\left(K{I}_1-\sum _{\alpha }{b}_{\alpha }{p}_{\alpha }-\alpha \omega \right)I+\left(2\mu -\frac{{\alpha }_J\omega }{J}\right)e^{\text{'}},$     (16)

${\phi }_1={\phi }_1^0+\left(b_1-{\phi }_1^0\right)I_1+\frac{p_1}{N_1}+\frac{p_2}{N_{12}}+{\alpha }_{p1}\omega ,$     (17)

${\phi }_2={\phi }_2^0+\left(b_2-{\phi }_2^0\right)I_1+\frac{p_2}{N_2}+\frac{p_1}{N_{12}}+{\alpha }_{p2}\omega ,$     (18)

где: $\lambda ,\mu$ – коэффициенты Ламе скелета; $K=\lambda +\frac{2}{3}\mu$ – модуль объемного сжатия скелета; $I_1=e:I,$ $J=\sqrt{e^{\text{'}}:e^{\text{'}}},$ $e^{\text{'}}$ – девиатор тензора малых деформаций; $b_{\alpha }$ и $N_{\alpha }$ – аналоги коэффициентов Био и модуля Био для матрицы и трещин; $N_{12}$ – модуль, отвечающий за взаимное влияние поровых давлений в подсистемах; $\alpha \ge \mathrm{0,}$ ${\alpha }_J\ge \mathrm{0,}$ ${\alpha }_{p1}\ge \mathrm{0,}$ ${\alpha }_{p2}\le 0$ – коэффициенты, отвечающие за разгрузку упругой энергии при развитии поврежденности; $\gamma ,$ $\beta$  – положительные коэффициенты, задающие нелинейную зависимость от поврежденности нулевого члена разложения упругой энергии скелета, характеризующего затраты энергии на появление новых поверхностей микротрещин в блоках матрицы.

Выражение для обменного члена в (1) и (2) в линейном приближении принимает вид:

$\dot{m}=\frac{{\rho }_{}^0}{{\mu }_f}A\left(\omega \right)\left(p_1-p_2\right).$     (19)

В случае мгновенной кинетики разрушения для параметра поврежденности справедливо выражение:

$\omega =\underset{\left[\mathrm{0,}t\right]}{\max }\left(\frac{1}{\beta }\left(\alpha I_1+{\alpha }_{J}J+{\alpha }_{p1}{p}_1+{\alpha }_{p2}{p}_2-\gamma \right)\right).$     (20)

Из выражения (20) следует критерий начала накопления поврежденности. Одновременное выполнение условий $\partial \omega /\partial t\ge 0$ и $\omega =0$ дает для границы зоны упругого поведения в пространстве инвариантов тензора малой деформации следующее уравнение:

$\alpha I_1+{\alpha }_{J}J+{\alpha }_{p1}{p}_1+{\alpha }_{p2}{p}_2-\gamma =0.$     (21)

Условие $\alpha I_1+{\alpha }_{J}J+{\alpha }_{p1}{p}_1+{\alpha }_{p2}{p}_2-\gamma <0$ соответствует упругому поведению. При выполнении условия $\alpha I_1+{\alpha }_{J}J+{\alpha }_{p1}{p}_1+{\alpha }_{p2}{p}_2-\gamma >0$ возможны активный процесс $\left(\partial \omega /\partial t>0\right)$ и упругая разгрузка (пассивный процесс). В пассивном процессе полагается $\omega =const.$

Подставляя определяющие соотношения (16)–(19) в (1)–(5), получим следующую систему уравнений:

$\frac{1}{M_1}\frac{\partial p_1}{\partial t}+\frac{1}{N_{12}}\frac{\partial p_2}{\partial t}+b_1\frac{\partial I_1}{\partial t}+{\alpha }_{p1}\frac{\partial \omega }{\partial t}=\frac{A\left(\omega \right)}{{\mu }_f}(p_2-p_1),$

(22)

$\begin{array}{c}\frac{1}{M_2}\frac{\partial p_2}{\partial t}+\frac{1}{N_{12}}\frac{\partial p_1}{\partial t}+b_2\frac{\partial I_1}{\partial t}+\\ +{\alpha }_{p2}\frac{\partial \omega }{\partial t}-\nabla \left(\frac{k_2\left(p_2\right)}{{\mu }_f}\nabla p_2\right)=\frac{A\left(\omega \right)}{{\mu }_f}(p_1-p_2),\end{array}$     (23)

$\left(\lambda +2\mu \right)\nabla I_1=b_1\nabla p_1+b_2\nabla p_2+\alpha \nabla \omega ,$    (24)

где $M_{\alpha }^{-1}=N_{\alpha }^{-1}+{\phi }_{\alpha }^{}{K}_f^{-1}.$

Волны депрессии и повреждения

Система уравнений (22)–(24) описывает сопряженные процессы движения флюидов, деформирования и рассеянного растрескивания скелета среды с двойной пористостью. Большой интерес представляет процесс растрескивания во время распространения волны депрессии в окрестности скважин или системы трещин гидроразрыва. Рассмотрим одномерный случай. Помимо относительной простоты этот случай имеет и практическое значение, поскольку он описывает процессы вблизи длинной плоской трещины гидроразрыва.

Действительно, основной метод добычи газа или нефти, разработанный для сланцевых месторождений в США, включает горизонтальное бурение и серию последовательных гидроразрывов пласта. В окрестности каждой трещины гидроразрыва создается область, называемая стимулированным объемом (stimulated reservoir volume, SRV), содержащая систему гидродинамически связанных вторичных трещин [Warpinski et al., 2009; Wu et al., 2014; Barati, Lang, 2014]. Протяженность стимулированного объема по нормали к плоскости трещин гидроразыва ограничена расстоянием между соседними трещинами, как показано на рис. 1, и может быть много меньше их длины (диаметра).

 

Рис. 1. К пояснению постановки задачи: 1 – скважина, 2 – магистральные трещины, 3 – стимулированный объем, содержащий систему вторичных трещин; X – протяженность стимулированного объема по нормали к плоскости трещины гидроразыва.

 

Рассмотрим плоскую трещину бесконечной протяженности в пористой среде с начальным давлением в каждой из подсистем $p_0.$ Пусть нормальная к плоскости трещины компонента полного тензора напряжения ${\sigma }_0,$ давление на плоской границе постоянно и меньше, чем $p_0.$ Пренебрежем силой гравитации. В этом случае уравнение равновесия (5) с учетом (16) после интегрирования по координате вдоль нормали к плоскости трещины принимает вид:

${\sigma }_0=\Lambda I_1-\sum _{\alpha }{b}_{\alpha }{p}_{\alpha }-\alpha \omega ,$     (25)

где $\Lambda =\lambda +2\mu .$ С учетом соотношения при одноосной деформации $J=-\sqrt{2/3}{I}_1$ и уравнения (25), соотношение (20) приводится к виду:

$\omega \left(x,t\right)=\underset{\tau \in \left[\mathrm{0,}t\right]}{\max }\left({\chi }_1{p}_1\left(x,\tau \right)+{\chi }_2{p}_2\left(x,\tau \right)-\Gamma ,0\right)$     (26)

где x – координата, в направлении нормали к плоскости трещины. В (26) использованы следующие обозначения:

$\Gamma =\frac{\gamma -\left(\alpha -{\alpha }_J\sqrt{2/3}\right){\sigma }_0/\Lambda }{\beta -\left(\alpha -{\alpha }_J\sqrt{2/3}\right)\alpha /\Lambda },$     (27)

$\begin{array}{c}{\chi }_1=\frac{{\alpha }_{p1}+\left(\alpha -{\alpha }_J\sqrt{2/3}\right)b_1/\Lambda }{\beta -\left(\alpha -{\alpha }_J\sqrt{2/3}\right)\alpha /\Lambda }\\ {\chi }_2=\frac{{\alpha }_{p2}+\left(\alpha -{\alpha }_J\sqrt{2/3}\right)b_2/\Lambda }{\beta -\left(\alpha -{\alpha }_J\sqrt{2/3}\right)\alpha /\Lambda }.\end{array}$     (28)

Согласно уравнению (26), если известны механические свойства материала, поврежденность зависит от истории изменения давления в подсистемах и напряженно-деформированного состояния среды. Далее рассматривается случай ${\chi }_2<\mathrm{0,}$ так как только в этом случае возможно развитие поврежденности. Пренебрегая зависимостью емкостных свойств матрицы от поврежденности, а также перекрестными членами в системе (22)–(23) по сравнению с диагональными членами, т.е., принимая условие:

$\frac{1}{N_{12}}+\frac{b_1{b}_2}{\Lambda }<<\frac{1}{M_k}+\frac{b_k^2}{\Lambda },k=\mathrm{1,}\mathrm{2,}$     (29)

получим следующую систему уравнений:

$\frac{\partial p_1}{\partial t}=A^{\text{'}}\left(\omega \right)\left(p_2-p_1\right),$     (30)

$\frac{\partial p_2}{\partial t}-\frac{\partial }{\partial x}\left(K_p\left(p_2\right)\frac{\partial p_2}{\partial x}\right)=r{A}^{\text{'}}\left(\omega \right)\left(p_1-p_2\right),$     (31)

где использованы обозначения:

$r=\left(\frac{1}{M_1}+\frac{b_1^2}{\Lambda }\right)/\left(\frac{1}{M_2}+\frac{b_2^2}{\Lambda }\right),$     (32)

$K_p=\frac{k_2\left(p_2\right)}{{\mu }_f}/\left(\frac{1}{M_2}+\frac{b_2^2}{\Lambda }\right),A^{\text{'}}\left(\omega \right)\frac{A\left(\omega \right)}{{\mu }_f}/\left(\frac{1}{M_1}+\frac{b_1^2}{\Lambda }\right).$

Система (30)–(31) решалась численно с использованием полностью неявной схемы второго порядка по пространству и первого порядка по времени.

Прежде обсуждения результатов численных расчетов, рассмотрим вариант задачи, допускающей аналитическое решение.

Автомодельное решение задачи о распространении волны депрессии и повреждения

Имеется связанный процесс депрессии и разрушения, вызванный понижением порового давления на границе полупространства x > 0 от начального пластового давления p₂(x, 0) = p₀, которое полагается равным в трещинах и матричных блоках, до некоторого значения p₂(0, t) = p_b. Для простоты считаем, что проводимость и емкость системы магистральных трещин – постоянные величины, зона трещиноватости является неограниченной. В начальный момент времени p₁(x, 0) = p₀, ω(x, 0) = 0. Процесс разрушения инициируется уменьшением давления p₂ до значения p_f, соответствующего началу разрушения в соответствии с (26). Выполнение условия $\omega =0$ определяет положение фронта разрушения. Поровое давление в матричных блоках до момента начала разрушения остается постоянным и равным начальному значению p0. После начала разрушения, вследствие появившейся гидродинамической связи между подсистемами двойной пористости, происходит постепенное выравнивание давлений p₁ и p₂. При $t\to \infty$ по мере того как протяженность ∆ пространственной области в которой происходит выравнивание давлений (∆ ~ $1/\sqrt{t})$ становится много меньше расстояния пройденного фронтом разрушения, решение для давлений как функция переменной  стремится к автомодельному пределу. В автомодельном решении зона выравнивания давлений представлена скачками давления p₁, производной $d{p}_2/d\xi$ и поврежденности на фронте разрушения $\xi ={\xi }_f.$ Уравнение для p₂ имеет вид:

$\frac{\partial p_2}{\partial t}-\frac{\partial }{\partial x}\left({\kappa }_2\frac{\partial p_2}{\partial x}\right)=\mathrm{0,}$ при $x>{\xi }_f\sqrt{t},$     (33)

$\left(1+r\right)\frac{\partial p_2}{\partial t}-\frac{\partial }{\partial x}\left({\kappa }^{*}\frac{\partial p_2}{\partial x}\right)=\mathrm{0,}$ при $x<{\xi }_f\sqrt{t}.$     (34)

Здесь $r={\phi }_1/{\phi }_2,$ ${\kappa }_2={\lambda }_2/\left({\phi }_{2}K\right),$ где: ${\lambda }_2,$ ${\phi }_2,$  – подвижность флюида и объемная доля (пористость) системы магистральных трещин; ${\kappa }^{*}={\lambda }^{*}/\left({\phi }_{2}K\right);$ ${\lambda }^{*}$ – подвижность флюида в разрушенной среде. Далее считаем, что ${\lambda }^{*}=const.$ Это возможно, если разрушение происходит только на фронте разрушения (это имеет место в случае ${\chi }_2<0$ и ${\chi }_1<0)$ или же, если (при ${\chi }_2<\mathrm{0,}$ ${\chi }_1>0)$  разрушение в области $x<{\Lambda }_f\sqrt{t}$ не приводит к существенному изменению подвижности флюида в системе магистральных трещин. В последнем случае можно считать ${\lambda }^{*}={\lambda }_2.$ Решение, удовлетворяющее условиям на разрыве, имеет вид:

$p_2(x,t)=\left\{\begin{array}{c}{p}_b+\frac{p_f-p_b}{erf\left(\frac{{\xi }_f}{\sqrt{4{\kappa }^{*}}}\right)}erf\left(\frac{x}{\sqrt{4{\kappa }^{*}t}}\right),x\le {\xi }_f\sqrt{t}\\ p_0-\frac{p_0-p_f}{1-erf\left(\frac{{\xi }_f}{\sqrt{4{\kappa }_2}}\right)}\left(1-erf\left(\frac{x}{\sqrt{4{\kappa }_{2}t}}\right)\right),x>{\xi }_f\sqrt{t}\end{array}\right.$

                                                  (35)

где ${\xi }_f$ определяется как решение нелинейного уравнения, которое следует из соотношений на разрыве:

$\begin{array}{c}\sqrt{{\kappa }^{*}}\frac{p_f-p_b}{erf\left(\frac{{\xi }_f}{\sqrt{4{\xi }^{*}}}\right)}\exp \left(-\frac{{\xi }_f^2}{4{\kappa }^{*}}\right)=\\ =\sqrt{{\kappa }_2}\frac{p_0-p_f}{1-erf\left(\frac{{\xi }_f}{\sqrt{4{\kappa }_2}}\right)}\exp \left(-\frac{{\xi }_f^2}{4{\kappa }_2}\right)+r\left(p_0-p_f\right){\xi }_f.\end{array}$     (36)

Решения, соответствующие различным значениям r, показаны на рис. 2.

 

Рис. 2. Асимптотическое при t → ∞ решение задачи о распространении волны разрушения в полупространстве x > 0, вызванной понижением порового давления в подсистеме трещин на его границе при x = 0. Показано отношение давлений p₂/p₀ как функции автомодельной переменной для четырех случаев, отличающихся значением r.

 

Излом на кривых соответствует положению фронта разрушения. При увеличении емкости матричных блоков (т.е. при увеличении r) скорость распространения фронта разрушения уменьшается. В соответствии с (35), зависимостью дебита q, приходящегося на единицу площади, от времени совпадает с зависимостью для случая отсутствия разрушения, а именно, $q~1/\sqrt{t}$ как в бесконечном однородном пласте. При этом на больших временах отношение величины дебита при наличии разрушения q к величине дебита при его отсутствии q₀ есть:

$\frac{q\left(t\right)}{q_0\left(t\right)}=\frac{p_f-p_b}{p_0-p_b}\frac{1}{erf\left({\xi }_f/\sqrt{4{\kappa }^{*}}\right)}.$     (37)

В зависимости от параметров задачи это отношение может быть произвольно большим. В этом случае зависимость дебита от времени при малых t имеет плато или рост на временном интервале, величина которого определяется временем перераспределения флюидов между подсистемами.

Результаты численных расчетов

Вернемся к полной постановке задачи. Пусть для простоты $K_p=const,$ $A^{\text{'}}\left(\omega \right)=A_{\omega }{\omega }^2.$  Для заданной протяженности зоны проводящих трещин X (см. рис. 1) определим временной масштаб $t_m=X^2/K_p.$ Используя масштабы $t_m,$ X, $p_0$ для давления, перепишем систему (30), (32) в безразмерном виде:

$\frac{\partial p_1}{\partial t}=A{\omega }^2\left(p_2-p_1\right),$     (38)

$\frac{\partial p_2}{\partial t}-\frac{{\partial }^2{p}_2}{\partial x^2}=rA{\omega }^2\left(p_1-p_2\right),$     (39)

где безразмерный коэффициент $A=A_{\omega }/t_m.$

При ${\chi }_2<0$ получается выражение для поврежденности:

$\omega \left(x,t\right)=\underset{\tau \in \left[\mathrm{0,}t\right]}{\max }\left(\chi p_1\left(x,\tau \right)-p_2\left(x,\tau \right)-{\Gamma }^{\text{'}},0\right).$     (40)

Здесь и далее имеются ввиду безразмерные величины, $\chi ={\chi }_1/\left|{\chi }_2\right|,$ ${\Gamma }^{\text{'}}=\Gamma /\left(\left|{\chi }_2\right|p_0^{}\right).$  Пусть

${\Gamma }^{\text{'}}=\chi p_1\left(x,0\right)-p_2\left(x,0\right).$     (41)

С учетом результатов (см. работу [Извеков и др., 2020]), справедливо приближенное равенство:

$\chi \approx \frac{1+\frac{K}{\Lambda }\left(1-\frac{\alpha }{{\alpha }_J}\sqrt{2/3}\right)b_1}{b_2|1-\frac{K}{\Lambda }\left(1-\frac{\alpha }{{\alpha }_J}\sqrt{2/3}\right)|}sign\left(\beta -\frac{{\alpha }^2}{\Lambda }\left(1-\frac{\alpha }{{\alpha }_J}\sqrt{2/3}\right)\right).$     (42)

Для сланцев можно сделать следующие оценки параметров модели [Извеков и др., 2020]. ${\alpha }_J/\alpha \approx 3-7\nu +2{\nu }^2/\sqrt{6}\left(1-{\nu }^2\right),$ $K/\Lambda =$$=3\left(1\text{\hspace{0.05em}}-2\nu \right)/\left(1\text{\hspace{0.05em}}-{\nu }^2\right),$ где $\nu$ – коэффициент Пуассона матрицы. Когда $\nu =0.2-0.35$ значение ${\alpha }_J/\alpha$ лежит в интервале 0.77–1.07, а $K/\Lambda \approx 1.88-\mathrm{1.03.}$ Величина $b_1/b_2\approx {10}^{-2}-{10}^2$ зависит от упругих свойств матрицы и характеристик системы трещин. Для определенности возьмем $b_1/b_2=\mathrm{10,}$ $b_1=\mathrm{1,}$ тогда параметр $\chi >0$ лежит в интервале 16.2–18.0 в зависимости от коэффициента Пуассона. Свойства горных пород варьируются в широком диапазоне, таким образом, нельзя утверждать, что параметр $\chi$ всегда положителен. Более того, система уравнений (38), (39) при положительных и при отрицательных $\chi$ имеет принципиально различные решения.

Дополним уравнения (38), (39) начальными условиями $p_1=p_2=\mathrm{1,}$ $\omega =0$  в слое $0<x<\mathrm{1,}$ и граничными условиями $p_2=0.1$ при $x=\mathrm{0,}$ $\partial p_2/\partial t=0$ при $x=1.$ Эта задача моделирует распространение волны депрессии и соответствующей зоны поврежденности в ограниченной зоне проводящих трещин в окрестности протяженной планарной трещины гидроразрыва. Решение рассматривается на отрезке времени, на котором волна депрессии не достигает правой границы слоя, т.е. изменением давления в проводящей подсистеме на правой границе слоя (x = 1) можно пренебречь. Примем в расчетах $A={10}^4.$ В зависимости от знака  возможны два типа решений. Случай $\chi <0.$ Повреждения начинают накапливаться в результате снижения давления в подсистеме 2. Снижение давления в подсистеме 1 ускоряет процесс разрушения. Эволюция зоны разрушения показана на рис. 3а. Случай $\chi >0.$ Как и в предыдущем случае, повреждения начинают накапливаться в результате снижения давления в подсистеме 2. Снижение давления в подсистеме 1 замедляет процесс разрушения. Результатом этого замедления является быстрое по сравнению с предыдущим случаем уменьшение скачка поврежденности в зависимости от x. Эволюция поврежденности показана на рис. 4а. Давления в подсистемах p₁ и p₂ в зависимости от расстояния до границы в момент времени t = 32 для рассматриваемых случаев показаны на рис. 3б и рис. 4б. Особенностью решения при $\chi >0$ является то, что, начиная с некоторого расстояния, за скачком поврежденности реализуется близкое к стационарному значение ω, что следует из сравнения значений параметра поврежденности при x = 0.2 в моменты времени t = 32 и t = 64.

 

Рис. 3. (a) – профиль поврежденности при $t=32$ (штрихпунктирная линия), при $t=64$ (сплошная линия); (б) – профили давлений при $t=32.$ Параметры ${\chi }_1=-\mathrm{0.5,}$ ${\chi }_2=-\mathrm{0.5,}$ $r={10}^3.$ 

 

Зависимость дебита от времени в широком диапазоне параметра r показана на рис. 4. Пунктирная линия соответствует решению без накопления повреждений, то есть без обмена массой между подсистемами. Из рисунка следует, что все решения имеют промежуточную асимптотику $q\propto 1/\sqrt{t}.$ Начальная стадия при больших r характеризуется практически постоянным во времени расходом. На больших временах видно снижение дебита из-за влияния ограниченной протяженности зоны проводящих трещин. При выбранном значении обменного коэффициента $A={10}^4$ при появлении поврежденности выравнивание давлений происходит достаточно быстро, таким образом, в рассматриваемом случае дебит не чувствителен к величине и знаку параметра $\chi .$

 

Рис. 4. (a) – профиль поврежденности при $t=32$  (штрихпунктирная линия), при $t=64$ (сплошная линия); (б) – профили давлений при $t=32.$  Параметры ${\chi }_1=\mathrm{0.5,}$ ${\chi }_2=-\mathrm{1.5,}$ $r={10}^3.$  

 

Рис. 5. Зависимость дебита от безразмерного времени для различных значений параметра r. Пунктирной линией показан дебит для неповреждаемой матрицы (без массообмена между подсистемами).

 

Обсуждение и выводы

Линеаризованная феноменологическая модель, описывающая сопряженные процессы течения жидкости в системе трещин, массообмен между матрицей и трещинами, деформацию скелета и повреждение матрицы содержит большое количество эмпирических параметров. В частной одномерной задаче модель сводится к более простой модели с меньшим количеством параметров (пороупругие константы и другие параметры объединяются в комплексы). Это позволяет прогнозировать распространение зоны поврежденности и влияние АВПД на дебит флюида. Постановка представленной в настоящей работе задачи приближенно соответствует конфигурации горизонтальной скважины с множественным гидроразрывом пласта. Предполагается, что в окрестности трещин гидроразыва предварительно созданы области гидродинамически связанной вторичной трещиноватости (стимулированный объем), протяженность которых ограничена расстоянием между магистральными трещинами. Задача при постоянной депрессии была решена численно. Решения имеют форму волны депрессии и следующей за ней волны разрушения. В начальный период дебит жидкости демонстрирует высокие значения, близкие к постоянным, что определяется параметрами модели. Все решения имеют промежуточную асимптотику, как в процессе без накопления повреждений. Такое поведение связано с предположением, что проницаемость системы трещин постоянна. Также получено автомодельное аналитическое решение в приближении скачкообразного накопления поврежденности на фронте разрушения.

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

Рассматриваемый процесс можно назвать “вынужденным авторазрывом”, так как причина разрыва – АВПД, а провоцируется процесс изменением внешних условий. Таким образом, собственно автофлюидоразрыв можно назвать “спонтанным авторазрывом”.

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

Финансирование работы

Исследование поддержано Российским научным фондом (грант №23-21-00175).

Об авторах

О.  Я.  Извеков

Московский физико-технический институт

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

А.  В.  Конюхов

Московский физико-технический институт; Объединенный институт высоких температур РАН

Email: izvekov_o@inbox.ru
Россия, Москва; Москва

Ю.  Н.  Извекова

Московский физико-технический институт

Email: izvekov_o@inbox.ru
Россия, Москва

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

Дополнительные файлы