Физика Земли, 2020, № 5, стр. 103-116
Термодинамически согласованная модель фильтрации в среде с двойной пористостью с учетом рассеянного разрушения матрицы
О. Я. Извеков 1, *, А. В. Конюхов 2, И. А. Чепрасов 1
1 Московский физико-технический институт (национальный исследовательский университет)
г. Москва, Россия
2 Объединенный институт высоких температур РАН
г. Москва, Россия
* E-mail: izvekov_o@inbox.ru
Поступила в редакцию 16.09.2019
После доработки 02.12.2019
Принята к публикации 27.01.2020
Аннотация
Развивается термодинамически согласованная модель среды с двойной пористостью с учетом рассеянного разрушения матрицы. Трещинообразованию в матрице способствует аномально высокое пластовое давление. Считается, что развитие поврежденности в матрице усиливает массообмен между подсистемами среды с двойной пористостью. На примере задачи о производительности длинной цилиндрической скважины качественно показано влияние аномально высокого пластового давления на дебит флюида.
ВВЕДЕНИЕ
Искусственные или естественные системы трещин играют ключевую роль в добыче нефти и газа из низкопроницаемых коллекторов [Olson, 2004]. Методы разработки таких месторождений, как правило, включают те или иные модификации гидравлического разрыва пласта (например, многостадийный гидроразрыв на горизонтальной скважине с использованием воды с добавлением присадок, снижающих ее вязкость) [Warpinski, 2009; Barati, 2014], в результате которых в окрестности скважины создается проницаемая область (SRV – stimulated reservoir volume) [Warpinski, 2009; Wu, 2014]. В настоящей работе изучается принципиальная возможность формирования подобной области в низкопроницаемых коллекторах, связанная только с аномально высоким пластовым давлением, и влияние ее формирования на дебит скважины.
Давление в порах называется аномально высоким, если оно превышает гидростатический уровень. Аномально высокое пластовое давление (АВПД) является распространенным явлением в пластовых породах со сверхнизкой проницаемостью, таких как карбонаты, известняки и сланцы [Fertl, 1976; Magara, 1978; Мелик-Пашаев, 1983; Luo, 2002]. Выделяют несколько причин возникновения АВПД, связанных с историей формирования залежи, термическим состоянием пласта, физико-химическими процессами, движением флюидов в земной коре. В частности, основным механизмом формирования АВПД в нефтематеринских породах является генерация флюида в условиях сверхнизкой проницаемости. АВПД достигает значений до 0.95 от горного давления [Li, 2012], и может не только служить причиной разрушения стенок скважин [Мелик-Пашаев, 1983; Li, 2012], но и способно поддерживать дебит нефти в низкопроницаемых коллекторах в отсутствие возможности применения традиционных методов поддержания пластового давления [Sonnenberg, 2009]. Кроме того, АВПД может приводить к развитию вторичных трещин в породе с низкой проницаемостью. Это явление, называемое автофлюидоразрыв (NHF – natural hydraulic fracturing), хорошо известно и обсуждается в связи с проблемой миграции флюидов [Secor, 1965; Engelder, 1990; Luo, 2002].
Модель двойной пористости [Barenblatt, 1960] широко используется при моделировании движения флюидов в трещиновато-пористых пластах [Голф-Рахт, 1986; Wu, 2014]. Трещиновато-пористые пласты состоят из блоков насыщенной пористой низкопроницаемой породы (матрицы), разделенных магистральными трещинами, по которым может двигаться жидкость. При этом матрица осуществляет подпитку системы трещин. В основе семейства континуальных моделей среды с двойной пористостью, насыщенной одной жидкостью, лежит гипотеза суперпозиции трех континуумов (предполагается, что элементарный объем, по которому проводится усреднение, много больше характерного расстояния между магистральными трещинами). Первый континуум образован частью жидкости, насыщающей матрицу, которая характеризуется низкой проницаемостью, но достаточно большой пористостью. Второй континуум образован жидкостью, заполняющей систему магистральных трещин, объемная доля которых мала по сравнению с пористостью матрицы. При этом проницаемость системы трещин предполагается достаточно высокой, чтобы обеспечить подвижность второго континуума. Третьим континуумом является твердый скелет, который может быть жестким, пороупругим или проявлять вязкие и пластические свойства. Между первым и вторым континуумом возможен обмен массой за счет перетекания жидкости из матрицы в трещины и наоборот. Если блоки матрицы не изолированы и проницаемость матрицы сравнима с проницаемостью трещин, модель двойной пористости часто называют моделью двойной проницаемости.
Как правило, учет влияния АВПД на производство нефти или газа в сланцевых породах сводится к учету зависимости проницаемости системы трещин в среде с двойной пористостью от давления [Thompson, 2010] в связи с тем, что для движения флюидов в низкопроницаемых породах требуются значительные депрессии. В отличие от указанного подхода в данной работе развивается модель среды с двойной пористостью, матрица которой способна растрескиваться под действием АВПД. Первоначально концепция автофлюидоразрыва была развита в работе [Secor, 1965] на основе подходов к моделированию традиционного гидроразрыва. В работе [Secor, 1965] рассматривается набор произвольно ориентированных трещин с непроницаемыми стенками на достаточно большом расстоянии друг от друга, чтобы исключить их взаимное влияние. В качестве критерия распространения трещин используется вариант энергетической теории Гриффитса. В дальнейшем эта модель была обобщена на случай пороупругой среды в предположении, что давление в трещинах в каждый момент времени равно поровому давлению в пространстве между трещинами [Engelder, 1990]. Указанные модели изначально сформулированы для описания развития трещиноватости в горных породах на геологических масштабах времени. При этом указывается, что автофлюидоразрыв играет важную роль в миграции флюидов. Например, если в силу определенных причин поровое давление в низкопроницаемом пласте, постепенно повышаясь, достигает критерия страгивания трещин (initiation), трещины начинают расти (propagation), что приводит к увеличению их объема, падению давления и, в результате, прекращению роста трещин (arrest). Процесс повторяется циклически. В работе [Luo, 2002] показано, что в сланцевых пластах при наличии АВПД возможно установление равновесия в окрестности критерия разрушения, когда генерируемый в нефтематеринской породе флюид отводится по системе трещин автофлюидоразрыва (зона I на рис. 1).
В задаче о производительности скважин характерные масштабы времени много меньше геологических, таким образом, возможной генерацией флюида в нефтематеринской породе можно пренебречь. В настоящей работе рассматривается ситуация, когда за длительный предшествующий период давление в матрице превысило гидростатический уровень, однако критерий разрушения еще не достигнут (например, зона II на рис. 1). Считается, что разрушение в матрице может произойти при участии АВПД в процессе технологических операций (например, бурение скважины, депрессия на забое скважины в процессе добычи и др.).
Развитие системы трещин, длина которых много меньше характерного размера элементарного объема, естественно моделировать методами теории повреждаемых сред (CDM – continuum damage mechanics) [Lemaitre, 1996; Murakami, 2012], в рамках которой коллектив микротрещин описывается осредненно с помощью специального параметра поврежденности. В настоящей работе для моделирования автофлюидоразрыва в матрице используется модель повреждаемости, аналогичная модели из работ [Kondaurov, 2009; Извеков, 2009; 2010]. Для простоты использован скалярный параметр поврежденности, что соответствует случайному распределению трещин по направлениям в пространстве. Критерий инициации начала развития повреждаемости в модели работ [Kondaurov, 2009; Извеков, 2009; 2010] является своеобразным обобщением энергетического критерия Гриффитса на случай рассеянного разрушения хрупких сред. Предполагается, что возникающие микротрещины в матрице в первую очередь усиливают массообмен между матрицей и магистральными трещинами.
ФОРМУЛИРОВКА МОДЕЛИ
Рассмотрим насыщенную пористую среду с упругим хрупким скелетом, включающим пористую матрицу и систему магистральных трещин. Матрица имеет низкую проницаемость и относительно высокую пористость. Система трещин характеризуется высокой проницаемостью и низкой пористостью (относительный объем пустот, содержащихся в трещинах). Матрица и трещины насыщены одной и той же слабосжимаемой жидкостью. Блоки матрицы изолированы друг от друга в том смысле, что непосредственное перемещение жидкости из блока в блок отсутствует и может осуществляться только через систему трещин.
Будем рассматривать среду с двойной пористостью как суперпозицию трех сплошных сред – двух флюидов и скелета. Первый континуум (флюид) соответствует жидкости, находящейся в матрице, второй – жидкости, заполняющей магистральные трещины. Жидкость взаимодействует со скелетом посредством массовых сил вязкого трения. В настоящей работе всюду используется изотермическое приближение (это предположение является традиционным при моделировании нетепловых методов добычи нефти). Предполагается наличие умеренного массообмена между матрицей и трещинами (это означает, что поток между матрицей и трещинами мал по сравнению с потоком в трещинах). Все движения предполагаются квазистатическими (что позволяет пренебречь кинетической энергией и влиянием массообмена на импульс), а деформации скелета малыми.
Далее используются индексы: “1”, “2”, “S” для обозначения величин, относящихся к флюиду, соответствующему жидкости в матрице, флюиду, соответствующему жидкости в трещинах, и к скелету соответственно. Также для обозначения величин, относящихся только к флюидам, используется индекс “α”, а к любому из континуумов индекс “A”. Определим величину ${{{\varphi }}_{{\alpha }}} = \frac{{dV_{{void}}^{{\alpha }}}}{{dV}}$ – объемную долю флюидов, где $dV_{{void}}^{{\alpha }}$ – объем, занимаемый жидкостью в подсистеме α (объем соответствующего пустотного пространства), $dV$ – элементарный объем. Определим плотности флюидов ${{r}_{{\alpha }}} = \frac{{d{{m}_{{\alpha }}}}}{{dV}}$ скелета ${{r}_{S}} = \frac{{d{{m}_{S}}}}{{dV}}$, где ${{m}_{{\alpha }}}$ ${{m}_{S}}$ – масса флюида в подсистеме α и масса скелета, содержащиеся в элементарном объеме $dV$, соответственно. Используя определение ${{\phi }_{{\alpha }}}$, получим, что ${{r}_{{\alpha }}} = {{\phi }_{{\alpha }}}{{{\rho }}_{f}}$ и ${{r}_{S}} = \left( {1 - {{\phi }_{1}} - {{\phi }_{2}}} \right){{{\rho }}_{S}}$, где ${{{\rho }}_{f}}$ – истинная плотность жидкости и ${{{\rho }}_{S}}$ – плотность материала скелета.
Выпишем основные уравнения в рамках сформулированных гипотез. Следствием закона сохранения массы являются три уравнения неразрывности, имеющие вид:
(1)
$\frac{{\partial {{r}_{1}}}}{{\partial t}} + \nabla \cdot \left( {{{r}_{1}}{{{\mathbf{v}}}_{1}}} \right) = {{\dot {m}}_{{12}}},$(2)
$\frac{{\partial {{r}_{2}}}}{{\partial t}} + \nabla \cdot \left( {{{r}_{2}}{{{\mathbf{v}}}_{2}}} \right) = {{\dot {m}}_{{21}}},$(3)
$\frac{{\partial {{r}_{S}}}}{{\partial t}} + \nabla \cdot \left( {{{r}_{S}}{{{\mathbf{v}}}_{S}}} \right) = 0,$Возьмем уравнения движения флюидов в форме закона Дарси
(4)
${{{\mathbf{W}}}_{\alpha }} = - \frac{{{{k}_{\alpha }}}}{{{{\mu }_{f}}}}\left( {\nabla {{p}_{\alpha }} - {{\rho }_{\alpha }}{\mathbf{g}}} \right),$Запишем уравнение равновесия для скелета и флюидов как единое целое:
Определим тензор Коши полных напряжений в пористой среде:
(6)
${\mathbf{T}} = {{{\mathbf{T}}}_{s}} - \left( {{{p}_{1}}{{\phi }_{1}} + {{p}_{2}}{{\phi }_{2}}} \right){\mathbf{I}},$Массовые силы взаимодействия между флюидами и скелетом, отнесенные к единице объема, имеют вид [Кондауров, 2007; Coussy, 2001]
(7)
${\mathbf{b}}_{\alpha }^{{\operatorname{int} }} = - {{p}_{\alpha }}\nabla {{\phi }_{\alpha }} + {\mathbf{b}}_{\alpha }^{{dis}},$Дифференциальный закон сохранения энергии для элементарного тела, состоящего из частиц трех континуумов, с учетом уравнений баланса массы (1)–(3) и уравнений (4), (5) имеет вид:
(8)
$\begin{gathered} \sum\limits_A {{{r}_{A}}\frac{{{{d}_{A}}{{u}_{A}}}}{{dt}}} = \sum\limits_A {{{{\mathbf{T}}}_{A}}:\nabla \otimes {{{\mathbf{v}}}_{A}}} - \sum\limits_A {{\mathbf{b}}_{\alpha }^{{\operatorname{int} }}} \cdot {{{\mathbf{w}}}_{\alpha }} + \\ + \,\,{{r}_{\Sigma }}Q - \nabla \cdot {\mathbf{q}} - \sum\limits_{\alpha \ne \beta } {{{{\dot {m}}}_{{\alpha \beta }}}{{u}_{\alpha }}} , \\ \end{gathered} $Для того же элементарного тела второе начало термодинамики в форме неравенства Клаузиуса–Дюгема с учетом массообмена имеет вид:
(9)
$\sum\limits_A {{{r}_{A}}\frac{{{{d}_{A}}{{\eta }_{A}}}}{{dt}} + } \sum\limits_{\alpha \ne \beta } {{{\eta }_{\alpha }}{{{\dot {m}}}_{{\alpha \beta }}}} + \nabla \left( {\frac{{\mathbf{q}}}{\theta }} \right) - \frac{{{{r}_{\Sigma }}Q}}{\theta } \geqslant 0,$Вводя удельную свободную энергию континуумов ${{\psi }_{A}} = {{u}_{A}} - {{\eta }_{A}}\theta $ и исключая из (9) поток приведенного тепла с помощью (8), а также используя (1)–(5), получим приведенное неравенство Клаузиуса–Дюгема в изотермическом приближении:
(10)
$\begin{gathered} - \sum\limits_A {{{r}_{A}}\frac{{{{d}_{A}}{{\psi }_{A}}}}{{dt}}} + {{{\mathbf{T}}}_{s}}:\frac{{{{d}_{s}}{\mathbf{e}}}}{{dt}} + \\ + \,\,\sum\limits_\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}} \geqslant 0, \\ \end{gathered} $Для замыкания системы законов сохранения (1)–(5) сформулируем систему термодинамически согласованных определяющих соотношений. Согласно теории, развитой Колеманом, Ноллом и Трусделлом [Трусделл, 1975], определяющие соотношения должны удовлетворять ряду утверждений, среди которых принцип локального действия (изменение состояния в некоторой точке не может непосредственно повлиять на реакцию материала на конечном от нее расстоянии), принцип инвариантности от системы отсчета (вид определяющего соотношения должен быть одним и тем же в любых системах отсчета – инерциальных и неинерциальных), а также принцип термодинамической согласованности. Последний принцип утверждает, что второе начало термодинамики в форме неравенства Клаузиуса–Дюгема (9) или (10) накладывает ограничение на возможную реакцию среды в произвольном процессе.
Определим понятия “процесс” и “реакция” в рассматриваемом случае двойной пороупругой среды с хрупким склетом. Сначала определим состояние материальной частицы пористой среды как набор следующих параметров: тензор малой деформации скелета ${\mathbf{e}}$, истинные плотности жидкостей в матрице и трещинах ${{\rho }_{\alpha }}$, относительная скорость флюидов и скелета ${{{\mathbf{w}}}_{\alpha }}$, параметр поврежденности в матрице $\omega $. Параметр $\omega $ является внутренней переменной, для которой необходимо сформулировать кинетическое (эволюционное) уравнение. Будем считать, что $\omega $ накапливается необратимо, т.е. залечивания микротрещин не происходит. Под “процессом” будем понимать последовательность смены параметров состояния, а под “реакцией” – значения остальных неизвестных физических величин.
В общем виде, в рамках сделанных предположений, в качестве определяющих соотношений рассматриваемой среды примем набор функций, связывающих реакцию $\Upsilon $ с текущим состоянием пористой среды:
(11)
$\Upsilon = \Upsilon \left( {{\mathbf{e}},{{\rho }_{\alpha }},{{{\mathbf{w}}}_{\alpha }},\omega } \right),$(12)
$\frac{{{{d}_{S}}\omega }}{{dt}} = \Omega \left( {{\mathbf{e}},{{\rho }_{\alpha }},{{{\mathbf{w}}}_{\alpha }},\omega } \right),$Определяющие соотношения, необходимые и достаточные для выполнения (10), в произвольном процессе имеют следующий вид:
(13)
${{\psi }_{\alpha }} = {{\psi }_{\alpha }}\left( {{{\rho }_{\alpha }}} \right);\,\,\,\,{{p}_{\alpha }} = \rho _{\alpha }^{2}\frac{{\partial {{\psi }_{\alpha }}}}{{\partial {{\rho }_{\alpha }}}};$(14)
${\mathbf{T}} = {{r}_{s}}\frac{{\partial \Psi }}{{\partial {\mathbf{e}}}};\,\,\,\,{{\phi }_{\alpha }} = - {{r}_{s}}\frac{{\partial \Psi }}{{\partial {{p}_{\alpha }}}}.$В выражениях (14) использован потенциал скелета $\Psi \left( {{{{\mathbf{e}}}_{s}},{{p}_{\alpha }},\omega } \right) = F - \sum\nolimits_\alpha {\frac{{{{p}_{\alpha }}}}{{{{r}_{s}}}}{{\phi }_{\alpha }}} $ [Кондауров, 2007], где $F\left( {{{{\mathbf{e}}}_{s}},{{\phi }_{\alpha }},\omega } \right)$ = ${{\psi }_{s}}\left( {{{{\mathbf{e}}}_{s}},{{\phi }_{\alpha }}\left( {{{{\mathbf{e}}}_{s}},{{\rho }_{\alpha }},\omega } \right),\omega } \right)$.
С учетом соотношений (13) и (14) неравенство (10) принимает вид:
где ${{\delta }_{\omega }} = - \left( {\frac{{\partial \Psi }}{{\partial \omega }}} \right)\frac{{{{d}_{s}}\omega }}{{dt}} = - \left( {\frac{{\partial \Psi }}{{\partial \omega }}} \right)\Omega $ – диссипация, связанная с растрескиванием матрицы.Процедура доказательства соотношений (13)–(15) во многом аналогична используемой в работе [Кондауров, 2007] и здесь не приводится.
Далее примем гипотезу, что каждая диссипация в отдельности неотрицательна, что является достаточным условием для выполнения (15):
ЛИНЕЙНАЯ МОДЕЛЬ
Линеаризуем полученные определяющие соотношения в окрестности особого состояния $\left\{ {0,\,\rho _{\alpha }^{0},\,0,\,\omega } \right\}$, т.е. примем, что в этом состоянии скелет не деформирован, жидкости покоятся, давления флюидов $p_{1}^{0} = p_{2}^{0} = 0$. Давления во флюидах и напряжения в скелете в произвольных состояниях будем считать малыми по сравнению с характерными упругими модулями скелета ${{{{p}_{\alpha }}} \mathord{\left/ {\vphantom {{{{p}_{\alpha }}} K}} \right. \kern-0em} K} \ll 1$, ${{\left\| {{{{\mathbf{T}}}_{S}}} \right\|} \mathord{\left/ {\vphantom {{\left\| {{{{\mathbf{T}}}_{S}}} \right\|} K}} \right. \kern-0em} K} \ll 1$, где $K$ – модуль объемного сжатия скелета. Также рассмотрим частный случай, в котором начальные, полные напряжения в состоянии $\left\{ {0,\,\rho _{\alpha }^{0},\,0,\,0} \right\}$, равны нулю (в общем случае это не так (см. работу [Мухамедиев, 2018]).
Разложение свободной энергии в ряд по малому параметру ${{\delta {{\rho }_{\alpha }}} \mathord{\left/ {\vphantom {{\delta {{\rho }_{\alpha }}} {\rho _{\alpha }^{0}}}} \right. \kern-0em} {\rho _{\alpha }^{0}}} \ll 1$ до квадратичных членов в окрестности состояния $\left\{ {0,\,\rho _{\alpha }^{0},\,0,\,\omega } \right\}{\text{:}}$
(17)
${{\rho }_{\alpha }}{{\psi }_{\alpha }}\left( {{{\rho }_{\alpha }}} \right) = \rho _{\alpha }^{0}\psi _{\alpha }^{0} + p_{\alpha }^{0}\left( {\frac{{\delta {{\rho }_{\alpha }}}}{{\rho _{\alpha }^{0}}}} \right) + \frac{1}{2}{{K}_{f}}{{\left( {\frac{{\delta {{\rho }_{\alpha }}}}{{\rho _{\alpha }^{0}}}} \right)}^{2}},$(18)
$\frac{{\delta {{\rho }_{\alpha }}}}{{\rho _{\alpha }^{0}}} = \frac{{\delta {{p}_{\alpha }}}}{{{{K}_{f}}}}.$Вслед за авторами работ [Kondaurov, 2009; Извеков, 2009; 2010] положим, что влияние поврежденности на упругие модули пренебрежимо мало. При этом макроскопическим проявлением поврежденности становится остаточная деформация в поврежденном материале, что формально сближает данный вариант модели повреждаемости с теорией пластичности. Скелет будем считать изотропным. Представим разложение потенциала скелета до величин второго порядка малости, считая поврежденность одним из малых параметров, в следующем виде:
(19)
$\begin{gathered} r_{s}^{0}\Psi ({\mathbf{e}},{{p}_{\alpha }},\omega ) = \frac{1}{2}\lambda I_{1}^{2} + \mu {{J}^{2}} + \sum\limits_\alpha {{{b}_{\alpha }}{{I}_{1}}{{p}_{\alpha }}} + \\ + \,\,\sum\limits_\alpha {\frac{1}{{2{{N}_{\alpha }}}}p_{\alpha }^{2}} + \frac{1}{{{{N}_{{12}}}}}{{p}_{1}}{{p}_{2}} + \sum\limits_\alpha {\phi _{\alpha }^{0}{{p}_{\alpha }}} + \gamma \omega + \\ + \,\,\beta {{{{\omega }^{2}}} \mathord{\left/ {\vphantom {{{{\omega }^{2}}} 2}} \right. \kern-0em} 2} - \alpha {{I}_{1}}\omega - {{\alpha }_{J}}J\omega - {{\alpha }_{{p1}}}{{p}_{1}}\omega - {{\alpha }_{{p2}}}{{p}_{2}}\omega , \\ \end{gathered} $Подстановка (19) в (14) дает линейные определяющие соотношения для скелета:
(20)
${\mathbf{T}} = \left( {K{{I}_{1}} - \sum\limits_\alpha {{{b}_{\alpha }}{{p}_{\alpha }} - \alpha \omega } } \right){\mathbf{I}} + \left( {2\mu - \frac{{{{\alpha }_{J}}\omega }}{J}} \right){\mathbf{e}}{\kern 1pt} ',$(21)
${{\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 {\kern 1pt} {\kern 1pt} ,$(22)
${{\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 {\kern 1pt} {\kern 1pt} ,$Обсудим неравенства диссипации (16). Учтем, что, согласно сделанным предположениям относительно начального состояния флюидов, $\psi _{1}^{0} = \psi _{2}^{0} = 0,$ тогда справедлива оценка:
(23)
$\begin{gathered} {{\delta }_{m}} = - \sum\limits_{\alpha \ne \beta } {\left( {{{\psi }_{\alpha }} + \frac{{{{p}_{\alpha }}}}{{{{\rho }_{\alpha }}}}} \right){{{\dot {m}}}_{{\alpha \beta }}}} \approx \\ \approx - \sum\limits_{\alpha \ne \beta } {\frac{{{{p}_{\alpha }}}}{{{{\rho }_{\alpha }}}}{{{\dot {m}}}_{{\alpha \beta }}}} \approx \frac{{\dot {m}}}{{\rho _{2}^{0}}}\left( {{{p}_{1}} - {{p}_{2}}} \right) \geqslant 0, \\ \end{gathered} $(24)
$\dot {m} = \frac{{\rho _{{}}^{0}}}{{{{\mu }_{f}}}}A\left( \omega \right)\left( {{{p}_{1}} - {{p}_{2}}} \right),$Поступая аналогично авторам [Kondaurov, 2009; Извеков, 2009; 2010], примем следующий закон эволюции поврежденности, достаточный для удовлетворения неравенству диссипации ${{\delta }_{\omega }} \geqslant 0$:
(25)
$\frac{{\partial \omega }}{{\partial t}} = - \frac{{{{r}_{S}}}}{{\tau \beta }}\left\langle {\frac{{\partial \Psi }}{{\partial \omega }}} \right\rangle ,\,\,\,\,\left\langle x \right\rangle = \left\{ \begin{gathered} x,\,\,x \geqslant 0, \hfill \\ 0,\,\,x < 0. \hfill \\ \end{gathered} \right.$Физический смысл выражения (25) заключается в том, что поврежденность развивается только в том случае, когда упругая энергия, высвобождающаяся при появлении новых поверхностей, превышает энергозатраты на их появление. Чем сильнее эта разница, тем интенсивнее развивается поврежденность. Уравнение (25) является в некотором роде аналогом теории Гриффитса для изолированной трещины.
Подставляя в (25) представление потенциала скелета (19), получим линейное эволюционное уравнение:
(26)
$\frac{{\partial \omega }}{{\partial t}} = \frac{1}{{\tau \beta }}\left\langle {\alpha {{I}_{1}} + {{\alpha }_{J}}J + {{\alpha }_{{p1}}}{{p}_{1}} + {{\alpha }_{{p2}}}{{p}_{2}} - \beta \omega - \gamma } \right\rangle {\text{.}}$Эволюция параметров однородного гидростатического напряженного состояния при отсутствии обмена между подсистемами описывается системой уравнений, следующих из (1), (2), (5) с учетом (20)–(22) и (26):
(28)
$\begin{gathered} \left( {\frac{1}{{{{N}_{1}}}} + \frac{{\phi _{1}^{0}}}{{{{K}_{F}}}} + \frac{{b_{1}^{2}}}{K}} \right)\frac{{\partial {{p}_{1}}}}{{\partial t}} + \left( {\frac{1}{{{{N}_{{12}}}}} + \frac{{{{b}_{1}}{{b}_{2}}}}{K}} \right)\frac{{\partial {{p}_{2}}}}{{\partial t}} + \\ + \,\,\left( {{{\alpha }_{{p1}}} + \frac{{\alpha {{b}_{1}}}}{K}} \right)\frac{{\partial \omega }}{{\partial t}} = 0, \\ \end{gathered} $(29)
$\begin{gathered} \left( {\frac{1}{{{{N}_{2}}}} + \frac{{\phi _{2}^{0}}}{{{{K}_{F}}}} + \frac{{b_{2}^{2}}}{K}} \right)\frac{{\partial {{p}_{2}}}}{{\partial t}} + \left( {\frac{1}{{{{N}_{{12}}}}} + \frac{{{{b}_{1}}{{b}_{2}}}}{K}} \right)\frac{{\partial {{p}_{1}}}}{{\partial t}} + \\ + \,\,\left( {{{\alpha }_{{p2}}} + \frac{{\alpha {{b}_{2}}}}{K}} \right)\frac{{\partial \omega }}{{\partial t}} = 0, \\ \end{gathered} $(30)
$\begin{gathered} \tau \frac{{\partial \omega }}{{\partial t}} = \left( {{{K}^{{ - 1}}}\frac{{{{\alpha }^{2}}}}{\beta } - 1} \right)\omega + \frac{1}{\beta }\left( {{{\alpha }_{{p1}}} + \frac{{\alpha {{b}_{1}}}}{K}} \right){{p}_{1}} + \\ + \,\,\frac{1}{\beta }\left( {{{\alpha }_{{p2}}} + \frac{{\alpha {{b}_{2}}}}{K}} \right){{p}_{2}} + \frac{1}{\beta }\left( {\frac{{\alpha С}}{K} - \gamma } \right). \\ \end{gathered} $Решение данной системы имеет вид ${\mathbf{x}} = a{{e}^{{\Gamma t}}} + {\mathbf{b}}$, где ${\mathbf{x}} = {{({{I}_{1}},{{p}_{1}},{{p}_{2}},\omega )}^{T}}$, $a = {\text{const}},\;{\mathbf{b}} = {\text{const}}$.
Система неравенств:
(32)
${{\left( {\frac{1}{{{{N}_{{12}}}}} + \frac{{{{b}_{1}}{{b}_{2}}}}{K}} \right)}^{2}} < \left( {\frac{1}{{{{N}_{1}}}} + \frac{{\phi _{1}^{0}}}{{{{K}_{F}}}} + \frac{{b_{1}^{2}}}{K}} \right)\left( {\frac{1}{{{{N}_{2}}}} + \frac{{\phi _{2}^{0}}}{{{{K}_{F}}}} + \frac{{b_{2}^{2}}}{K}} \right)$Из ограниченности скорости ${{\partial \omega } \mathord{\left/ {\vphantom {{\partial \omega } {\partial t}}} \right. \kern-0em} {\partial t}}$ при $\tau \to 0$, из (26) следует выражение для параметра $\omega $:
(33)
$\omega = \mathop {\max }\limits_{\left[ {0,\,t} \right]} \left( {\frac{1}{\beta }\left( {\alpha {{I}_{1}} + {{\alpha }_{J}}J + {{\alpha }_{{p1}}}{{p}_{1}} + {{\alpha }_{{p2}}}{{p}_{2}} - \gamma } \right)} \right).$Из выражения (26) следует критерий начала накопления поврежденности. Одновременное выполнение условий ${{\partial \omega } \mathord{\left/ {\vphantom {{\partial \omega } {\partial t}}} \right. \kern-0em} {\partial t}} \geqslant 0$ и $\omega = 0$ дает для границы зоны упругого поведения в пространстве инвариантов тензора малой деформации следующее уравнение:
(34)
$\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$ соответствует упругому поведению. При выполнении условия $\alpha {{I}_{1}} + {{\alpha }_{J}}J + {{\alpha }_{{p1}}}{{p}_{1}} + {{\alpha }_{{p2}}}{{p}_{2}} - \gamma > 0$ возможны активный процесс (${{\partial \omega } \mathord{\left/ {\vphantom {{\partial \omega } {\partial t}}} \right. \kern-0em} {\partial t}} > 0$) и упругая разгрузка (пассивный процесс). В пассивном процессе полагается $\omega = {\text{const}}$.
ОЦЕНКА ПОРОУПРУГИХ КОЭФФИЦИЕНТОВ
Выяснению физического смысла и способам лабораторного определения пороупругих констант для среды с двойной пористостью посвящен ряд работ [Bai et al., 1995; Berryman, Pride, 2002]. Существуют также обобщения пороупругой модели на случай мультипористости [Mehrabian, Abousleiman, 2015]. Приведем здесь простые оценки относительных значений коэффициентов ${{b}_{1}}$ и ${{b}_{2}}$, модулей ${{N}_{\alpha }}$ и ${{N}_{{12}}}$, подвергнув среду одноосному сжатию с нулевой поперечной деформацией. Соответствующая приближенная одномерная эквивалентная схема состоит из двух последовательно соединенных герметичных камер (элементов схемы). В ненапряженном состоянии продольные размеры камер равны $l_{1}^{0}$ и $l_{2}^{0}$, причем $l_{2}^{0} \ll l_{1}^{0}$, как показано на рис. 2. Торцевые стенки камер могут скользить без трения вдоль горизонтальной оси. Первая камера, играющая роль матрицы, заполнена пороупругой средой с пористостью ${{n}_{1}}$, упругим модулем ${{E}_{1}}$, коэффициентом Био ${{\alpha }_{1}}$ и модулем Био $\tilde {N}$. Внутри второй камеры, соответствующей магистральным трещинам, находится пружина жесткостью ${{E}_{2}}$ (считаем, что трещины не заполнены ничем, кроме жидкости). Обозначим ${{\xi }_{\alpha }} = {{l_{\alpha }^{{}}} \mathord{\left/ {\vphantom {{l_{\alpha }^{{}}} {\left( {l_{1}^{0} + l_{2}^{0}} \right)}}} \right. \kern-0em} {\left( {l_{1}^{0} + l_{2}^{0}} \right)}}$ объемные доли матрицы и трещин, так что ${{\xi }_{1}} + {{\xi }_{2}} = 1$. Объемные доли флюидов (“пористости”) будут равны ${{\phi }_{1}} = {{n}_{1}}{{\xi }_{1}}$, ${{\phi }_{2}} = {{\xi }_{2}}$. В камеры независимо нагнетаются давления ${{p}_{1}}$ и ${{p}_{2}}$.
Запишем условия равновесия:
(35)
$\sigma = - {{\alpha }_{1}}{{p}_{1}} + {{E}_{1}}{{\varepsilon }_{1}} = - {{p}_{2}} + {{E}_{2}}{{\varepsilon }_{2}},$Запишем общую деформацию системы ${{\varepsilon }_{\Sigma }} = $ $ = {{\left( {\Delta {{l}_{1}} + \Delta {{l}_{2}}} \right)} \mathord{\left/ {\vphantom {{\left( {\Delta {{l}_{1}} + \Delta {{l}_{2}}} \right)} {\left( {l_{1}^{0} + l_{2}^{0}} \right)}}} \right. \kern-0em} {\left( {l_{1}^{0} + l_{2}^{0}} \right)}}$ как комбинацию относительных удлинений камер:
(36)
${{\varepsilon }_{\Sigma }} = {{\xi }_{1}}{{\varepsilon }_{1}} + {{\xi }_{2}}{{\varepsilon }_{2}}.$Выражая относительные удлинения из (36), с помощью (35) получим аналог закона Гука:
(37)
$\sigma = {{E}_{\Sigma }}{{\varepsilon }_{\Sigma }} - \left( {{{\alpha }_{1}}\frac{{{{E}_{\Sigma }}}}{{{{E}_{1}}}}{{\xi }_{1}}} \right){{p}_{1}} - \left( {\frac{{{{E}_{\Sigma }}}}{{{{E}_{2}}}}{{\xi }_{2}}} \right){{p}_{2}},$где $E_{\Sigma }^{{}} = {{\left( {{{{{\xi }_{1}}} \mathord{\left/ {\vphantom {{{{\xi }_{1}}} {{{E}_{1}}}}} \right. \kern-0em} {{{E}_{1}}}} + {{{{\xi }_{2}}} \mathord{\left/ {\vphantom {{{{\xi }_{2}}} {{{E}_{2}}}}} \right. \kern-0em} {{{E}_{2}}}}} \right)}^{{ - 1}}}$ – суммарный упругий модуль системы.
Сравнивая (37) с (20) в отсутствие поврежденности, приходим к выводу, что коэффициенты при давлениях ${{p}_{1}}$ и ${{p}_{2}}$ играют роль коэффициентов ${{b}_{1}}$ и ${{b}_{2}}$. Таким образом,
(38)
$\frac{{{{b}_{1}}}}{{{{b}_{2}}}} = \frac{{{{\alpha }_{1}}{{\xi }_{1}}}}{{{{\xi }_{2}}}}\frac{{{{E}_{2}}}}{{{{E}_{1}}}}.$При стремлении объемной доли ${{\xi }_{2}}$ к нулю, суммарный упругий модуль системы ${{E}_{\Sigma }}$ стремится к ${{E}_{1}}$, а коэффициент ${{b}_{1}}$ стремится к коэффициенту Био матрицы ${{\alpha }_{1}}$. Пусть ${{\xi }_{2}} \approx 0.01{{\xi }_{1}}$, тогда при соотношении модулей ${{{{E}_{2}}} \mathord{\left/ {\vphantom {{{{E}_{2}}} {{{E}_{1}}}}} \right. \kern-0em} {{{E}_{1}}}} \approx {{10}^{{ - 4}}} - 1$, отношение ${{{{b}_{1}}} \mathord{\left/ {\vphantom {{{{b}_{1}}} {{{b}_{2}}}}} \right. \kern-0em} {{{b}_{2}}}}$ будет лежать в диапазоне ${{10}^{{ - 2}}}{\kern 1pt} - {\kern 1pt} {{10}^{2}}$.
Модуль ${{E}_{2}}$ определяется как ${{E}_{2}} = {{S}_{n}}\Delta \lambda = {{S}_{n}}{{\xi }_{2}}L$, где: ${{S}_{n}}$ – нормальная жесткость одиночной трещины; $\Delta \lambda $ – ширина одиночной трещины; $L$ – характерный размер блока [Кочарян, 2016]. Возьмем ${{S}_{n}} \approx {{10}^{2}}$ ГПа/м, $L \approx 1$ м, ${{E}_{1}} = 10$ ГПа, согласующиеся с данными из табл. 1, получим, что отношение упругих модулей будет равно ${{{{E}_{2}}} \mathord{\left/ {\vphantom {{{{E}_{2}}} {{{E}_{1}}}}} \right. \kern-0em} {{{E}_{1}}}} \approx 0.1$, тогда согласно (38) отношение ${{{{b}_{1}}} \mathord{\left/ {\vphantom {{{{b}_{1}}} {{{b}_{2}}}}} \right. \kern-0em} {{{b}_{2}}}} \approx 10$, а из (37) следует ${{b}_{1}} = 0.9{{\alpha }_{1}}$, ${{b}_{2}} = 0.09{{\alpha }_{1}}$.
Таблица 1
Величина | Значение | Материал | Литературный источник |
---|---|---|---|
${{K}_{f}}$, ГПа | 0.5–5 | Нефть | [Craft, Hawkins, 1959] |
1.7 | [Abousleiman, Nguyen, 2005] | ||
${{K}_{1}}$, ГПа | 1.1 | Сланцы (Gulf of Mexico) | [Abousleiman, Nguyen, 2005] |
9.8 | Сланцы (Woodford shale) | [Abousleiman et al., 2007] | |
${{\alpha }_{1}}$ | 0.96 | Сланцы (Gulf of Mexico) | [Abousleiman, Nguyen, 2005] |
${{\xi }_{2}}$ | 0.0001–0.01 | [Snow, 1968] | |
${{S}_{n}}$, ГПа/м | 105–117 | Сланцы (Barnett) | [Ye et al., 2016] |
226 | Сланцы (Mancos) | [Ye et al., 2016] | |
25–40 | Сланцы (Pierre) | [Ye et al., 2016] | |
$L$, м | 0.2–30 | [Snow, 1968] | |
${{E}_{1}} = {{K}_{1}} + \frac{4}{3}\mu $, ГПа | 2.1 | Сланцы (Gulf of Mexico) | [Abousleiman, Nguyen, 2005] |
10–20 (анизотропия) | Сланцы (Woodford shale) | [Abousleiman et al., 2007] | |
5.1 ± 0.5 | Сланцы с высоким содержанием органики | [Eseme et al., 2007] | |
16.6 ± 2 | Сланцы с низким содержанием органики | [Eseme et al., 2007] | |
${{n}_{1}}$ | 3.8% | Сланцы (Woodford shale) | [Abousleiman et al., 2007] |
14% | Сланцы (Gulf of Mexico) | [Abousleiman, Nguyen, 2005] | |
$\tilde {N}$, ГПа | 27.57 | Сланцы (Gulf of Mexico) | [Abousleiman, Nguyen, 2005] |
$\nu $ | 0.11–0.29 (анизотропия) | Сланцы (Woodford shale) | [Abousleiman et al., 2007] |
0.35 | Сланцы с высоким содержанием органики | [Eseme et al., 2007] | |
0.2 | Сланцы с низким содержанием органики | [Eseme et al., 2007] |
Оценим модули ${{N}_{\alpha }}$ и ${{N}_{{12}}}$ с помощью предложенной упрощенной схемы. Учитывая, что деформация трещин ${{\varepsilon }_{2}} = {{\Delta {{\phi }_{2}}} \mathord{\left/ {\vphantom {{\Delta {{\phi }_{2}}} {{{\phi }_{2}}}}} \right. \kern-0em} {{{\phi }_{2}}}}$, деформация матрицы ${{\varepsilon }_{1}} = {{\Delta {{\xi }_{1}}} \mathord{\left/ {\vphantom {{\Delta {{\xi }_{1}}} {{{\xi }_{1}}}}} \right. \kern-0em} {{{\xi }_{1}}}} = {{\Delta {{\phi }_{1}}} \mathord{\left/ {\vphantom {{\Delta {{\phi }_{1}}} {{{\phi }_{1}}}}} \right. \kern-0em} {{{\phi }_{1}}}} - {{\Delta {{n}_{1}}} \mathord{\left/ {\vphantom {{\Delta {{n}_{1}}} {{{n}_{1}}}}} \right. \kern-0em} {{{n}_{1}}}}$, при этом $\Delta {{n}_{1}} = {{\alpha }_{1}}{{\varepsilon }_{1}} + {{{{p}_{1}}} \mathord{\left/ {\vphantom {{{{p}_{1}}} {\tilde {N}}}} \right. \kern-0em} {\tilde {N}}}$ из уравнений (35) и (36) получим оценку:
(39)
$\frac{1}{{{{N}_{1}}}} = \frac{{\left( {{{\alpha }_{1}} - {{b}_{1}}} \right)\left( {1 + {{{{\alpha }_{1}}} \mathord{\left/ {\vphantom {{{{\alpha }_{1}}} {{{n}_{1}}}}} \right. \kern-0em} {{{n}_{1}}}}} \right){{\phi }_{1}}}}{{{{E}_{1}}}} + \frac{{{{\phi }_{1}}}}{{{{n}_{1}}\tilde {N}}},$Тогда для относительных величин справедливо:
(42)
$\frac{{{{N}_{2}}}}{{{{N}_{1}}}} = \frac{1}{{1 - {{b}_{2}}}}\left( {\left( {{{\alpha }_{1}} - {{b}_{1}}} \right)\left( {1 + \frac{{{{\alpha }_{1}}}}{{{{n}_{1}}}}} \right) + \frac{{{{E}_{1}}}}{{{{n}_{1}}\tilde {N}}}} \right)\frac{{{{b}_{1}}}}{{{{b}_{2}}}},$При ${{{{b}_{1}}} \mathord{\left/ {\vphantom {{{{b}_{1}}} {{{b}_{2}}}}} \right. \kern-0em} {{{b}_{2}}}} \approx 10$ и ${{\alpha }_{1}} = 0.7$ из (43) получим оценку ${{{{N}_{2}}} \mathord{\left/ {\vphantom {{{{N}_{2}}} {{{N}_{{12}}}}}} \right. \kern-0em} {{{N}_{{12}}}}} \approx 0.68$. Если дополнительно положить ${{n}_{1}} = 0.1$ и ${{{{E}_{1}}} \mathord{\left/ {\vphantom {{{{E}_{1}}} {\left( {{{n}_{1}}\tilde {N}} \right)}}} \right. \kern-0em} {\left( {{{n}_{1}}\tilde {N}} \right)}} \approx 1$, тогда согласно соотношению (42) отношение ${{{{N}_{2}}} \mathord{\left/ {\vphantom {{{{N}_{2}}} {{{N}_{1}}}}} \right. \kern-0em} {{{N}_{1}}}} \approx 16.6$. Таким образом, абсолютные значения модулей ${{N}_{2}}$ и ${{N}_{{12}}}$ в рамках сделанных предположений имеют один порядок и более чем на порядок превышают модуль ${{N}_{1}}$. В последующих расчетах будем придерживаться полученных соотношений.
ОЦЕНКА ПАРАМЕТРОВ В МОДЕЛИ ПОВРЕЖДЕННОСТИ ДЛЯ МАТРИЦЫ
В предложенную линейную модель поврежденности входит ряд коэффициентов. Оценим их, задавшись гипотезой о частном механизме разрушения матрицы. Обозначим коэффициенты, отвечающие за уменьшение упругой энергии матрицы ${{\alpha }_{M}}$, ${{\alpha }_{{JM}}}$, ${{\alpha }_{{pM}}}$. В работах [Engelder, 1990; Luo, 2002] указывается, что трещины автофлюидоразрыва развиваются преимущественно за счет нормального раскрытия. Принимая указанный механизм разрушения, положим, что в матрице содержится некоторое количество плоских круглых трещин (penny shaped) диаметром $l$, ориентация которых имеет изотропное распределение. Пренебрежем взаимодействием этих трещин, таким образом, уравнение (34) при ${{p}_{2}} = 0$ соответствует условию начала развития одиночной трещины. В отличие от работы [Engelder, 1990] будем полагать, что стенки трещин непроницаемы для флюида, что правомерно в условиях сверхнизкой проницаемости.
В случае одноосного сжатия (рис. 3а) под действием внутреннего давления первой раскроется трещина, ориентированная перпендикулярно максимальному сжимающему напряжению ${{\sigma }_{V}} < 0$. Это произойдет при внутреннем давлении ${{p}_{{ini}}}$, равном [Engelder, 1990]:
(44)
${{p}_{{ini}}} = \frac{{{{K}_{{IC}}}}}{{1.13\sqrt l }} - \frac{\nu }{{1 - \nu }}{{\sigma }_{V}},$В случае всестороннего сжатия (рис. 3б) уравнение для критического давления будет иметь вид:
Сравнивая (28) с (47), получим оценки:
Аналогично рассматривая трещину в условиях одноосного сжатия (в этом случае $J = - \sqrt {{2 \mathord{\left/ {\vphantom {2 3}} \right. \kern-0em} 3}} {{I}_{1}}$), получим:
(48)
$\frac{{{{\alpha }_{{JM}}}}}{{{{\alpha }_{M}}}} = \frac{{3 - 7\nu + 2{{\nu }^{2}}}}{{\sqrt 6 {{{\left( {1 - \nu } \right)}}^{2}}}}.$В работе [Grytsenko, Galybin, 2010] с помощью численных расчетов методом сингулярных интегральных уравнений было показано, что в случае ансамбля из большого числа одинаковых трещин, находящихся на расстояниях больше либо порядка их длины, критический коэффициент интенсивности напряжений зависит от взаимного расположения трещин, однако отличается от ${{K}_{{IC}}}$ не более чем на 30% (в меньшую сторону).
Параметр $\beta $ в расчетах подбирался таким, чтобы выполнялось неравенство ${{\alpha _{M}^{2}} \mathord{\left/ {\vphantom {{\alpha _{M}^{2}} \beta }} \right. \kern-0em} \beta } < K.$
УСЛОВИЕ РАЗРУШЕНИЯ СРЕДЫ С ДВОЙНОЙ ПОРИСТОСТЬЮ
Обратимся теперь к среде с двойной пористостью. Воспользуемся одномерной моделью. Пусть условие начала поврежденности в матрице имеет вид:
Требуется переформулировать его в терминах полной деформации ${{\varepsilon }_{\Sigma }}$. Исключим ${{\varepsilon }_{1}}$ из условия (49) с помощью (35) и (36). В итоге получим:
(50)
$\begin{gathered} \frac{{{{\alpha }_{M}}}}{{{{\xi }_{1}}}}\left( {1 - {{b}_{2}}} \right){{\varepsilon }_{\Sigma }} + {{\alpha }_{M}}\left( {{{a}_{1}}{{p}_{1}} - {{a}_{2}}{{p}_{2}}} \right) + \\ + \,\,{{\alpha }_{p}}_{M}{{p}_{1}} - \gamma = 0, \\ \end{gathered} $Сравнивая полученное выражение с (34), получим:
Используя ранее сделанные оценки ${{{{b}_{1}}} \mathord{\left/ {\vphantom {{{{b}_{1}}} {{{b}_{2}}}}} \right. \kern-0em} {{{b}_{2}}}} \approx 10$, ${{{{E}_{2}}} \mathord{\left/ {\vphantom {{{{E}_{2}}} {{{E}_{1}}}}} \right. \kern-0em} {{{E}_{1}}}} \approx 0.1$, ${{\xi }_{2}} \approx 0.01{{\xi }_{1}}$, получим ${{a}_{1}} \approx {{a}_{2}}$, $\alpha \approx {{\alpha }_{M}}$. Далее, используя оценку $\frac{{{{\alpha }_{M}}}}{{{{\alpha }_{p}}_{M}}} = {{E}_{1}}$ (одномерный аналог (46)), получим соотношение между коэффициентами ${{\alpha }_{{p1}}}$ и ${{\alpha }_{{p2}}}$:
ЗАДАЧА О ПРОИЗВОДИТЕЛЬНОСТИ СКВАЖИНЫ
Далее для определенности пренебрежем силой тяжести. Рассматривается однородный безграничный пласт. Начальные давления в трещинах и в матрице пласта одинаковы и равны ${{p}_{0}}$. На бесконечности задано сжимающее гидростатическое горное давление $\sigma _{0}^{{}} < 0$. В некоторый момент времени в пласте появляется длинная цилиндрическая скважина радиусом $a$, давление на забое которой устанавливается на заданном уровне ниже начального давления. В результате напряженное состояние среды изменяется, а флюиды из системы трещин начинают поступать в скважину. В процессе развития зоны депрессии в окрестности скважины под действием АВПД матрица растрескивается, что усиливает массообмен с трещинами.
Решим задачу численно в приближении радиальной симметрии, в котором перемещение частиц скелета происходит только в радиальном направлении и не зависит от полярного угла. Обозначим радиальную компоненту перемещения $u$. Тогда тензор напряжений в скелете принимает вид:
(55)
${\mathbf{e}} = \frac{{\partial u}}{{\partial r}}{{{\mathbf{e}}}_{r}} \otimes {{{\mathbf{e}}}_{r}} + \frac{u}{r}{{{\mathbf{e}}}_{\theta }} \otimes {{{\mathbf{e}}}_{\theta }},$(57)
$J = \sqrt {{\mathbf{e}}{\kern 1pt} ':{\mathbf{e}}{\kern 1pt} '} = \sqrt {\frac{2}{3}} \sqrt {{{{\left( {\frac{{\partial u}}{{\partial r}} - \frac{u}{r}} \right)}}^{2}} + \frac{u}{r}\frac{{\partial u}}{{\partial r}}{\kern 1pt} } .$Подставим определяющие соотношения (20)–(22), закон массобмена (24) и закон Дарси (4) в уравнение равновесия (5) и уравнения неразрывности (1), (2). Получим следующую систему уравнений:
(58)
$\begin{gathered} \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(\omega )}}{{{{\mu }_{f}}}}({{p}_{2}} - {{p}_{1}}), \\ \end{gathered} $(59)
$\begin{gathered} \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}} - \\ - \,\,\frac{1}{r}\frac{\partial }{{\partial r}}\left( {\frac{{{{k}_{2}}({{p}_{2}})}}{{{{\mu }_{f}}}}r\frac{{\partial {{p}_{2}}}}{{\partial r}}} \right) = \frac{{A(\omega )}}{{{{\mu }_{f}}}}({{p}_{1}} - {{p}_{2}}), \\ \end{gathered} $(60)
$(\lambda + 2\mu )\frac{{\partial {{I}_{1}}}}{{\partial r}} = {{b}_{1}}\frac{{\partial {{p}_{1}}}}{{\partial r}} + {{b}_{2}}\frac{{\partial {{p}_{2}}}}{{\partial r}} + \alpha \frac{{\partial \omega }}{{\partial r}},$Для постановки граничных условий заменим бесконечный пласт расчетной областью цилиндрической формы, так чтобы оси скважины и расчетной области совпадали. Пусть радиус расчетной области $R$ много больше радиуса скважины. Изначально недеформированную расчетную область (без скважины) подвергнем гидростатическому сжимающему напряжению $\sigma _{0}^{{}} < 0$ при условии постоянства давлений $p_{0}^{{}} > 0$ в матрице и трещинах. Перемещение границы на расстоянии $R \gg a$ от скважины будет равно:
(61)
${{u}_{R}} = \frac{{{{\sigma }_{0}} + ({{b}_{1}} + {{b}_{2}}){{p}_{0}}}}{{2(\lambda + \mu )}}R \approx \frac{{{{\sigma }_{0}} + {{b}_{1}}{{p}_{0}}}}{{2(\lambda + \mu )}}R.$В расчетах будем считать эту границу закрепленной, т.е. $u(R,\,t) = {{u}_{R}}$. Поставим следующие граничные условия на давление в трещинах:
Условие равенства полных радиальных напряжений на стенках скважины $\left( {r = a} \right)$ имеет вид:
(64)
$ - {{p}_{a}} = \lambda \left[ {{{{\left( {\frac{{\partial u}}{{\partial r}}} \right)}}_{a}} + \frac{{{{u}_{a}}}}{a}} \right] + 2\mu {{\left( {\frac{{\partial u}}{{\partial r}}} \right)}_{a}} - {{b}_{1}}{{p}_{{1a}}} - {{b}_{2}}{{p}_{a}}.$В качестве начальных условий положим условие на начальное поле давлений
(65)
${{p}_{1}}\left( {r,\,0} \right)\,\,{\text{ = }}\,\,{{p}_{2}}\left( {r,\,0} \right) = {{p}_{0}}$(66)
$u\left( {r,\,0} \right) = \frac{{{{\sigma }_{0}} + ({{b}_{1}} + {{b}_{2}}){{p}_{0}}}}{{2(\lambda + \mu )}}r,\,\,\,\,r \geqslant a$СИСТЕМА УРАВНЕНИЙ В БЕЗРАЗМЕРНОМ ВИДЕ
Для представления системы (58)–(60) в безразмерном виде введем характерные масштабы для давления, напряжения и упругих модулей $\mu $, для координат и перемещений ${{r}_{0}}$, для времени ${{t}_{0}}$, для проницаемости системы трещин $k_{2}^{0}$. Увяжем масштабы ${{t}_{0}}$ и $k_{2}^{0}$ с помощью следующего соотношения:
Безразмерную функцию поврежденности, отвечающую за обмен между подсистемами трещины и матрицы, представим в виде $\overline {A(\omega )} = {{r_{0}^{2}} \mathord{\left/ {\vphantom {{r_{0}^{2}} {k_{2}^{0}}}} \right. \kern-0em} {k_{2}^{0}}}A(\omega )$. Безразмерные коэффициенты в выражении для поврежденности $\bar {\alpha } = {\alpha \mathord{\left/ {\vphantom {\alpha \mu }} \right. \kern-0em} \mu },$ ${{\bar {\alpha }}_{p}} = {{\alpha }_{p}},$ $\bar {\beta } = {\beta \mathord{\left/ {\vphantom {\beta \mu }} \right. \kern-0em} \mu },$ $\bar {\gamma } = {\gamma \mathord{\left/ {\vphantom {\gamma \mu }} \right. \kern-0em} \mu }.$ Безразмерные коэффициенты Ламе $\bar {\lambda } = {\lambda \mathord{\left/ {\vphantom {\lambda \mu }} \right. \kern-0em} \mu }$, $\bar {\mu } = 1$.
После обезразмеривания система (58)–(60) принимает следующий вид (все переменные далее безразмерные, черточки над ними для удобства опущены):
(68)
$\begin{gathered} \frac{{\partial {{p}_{1}}}}{{\partial t}} + {{\varepsilon }_{M}}{{\varepsilon }_{{12}}}\frac{{\partial {{p}_{2}}}}{{\partial t}} + {{b}_{1}}\varepsilon {{\varepsilon }_{M}}\frac{\partial }{{\partial t}}\left( {\frac{{\partial u}}{{\partial r}} + \frac{u}{r}} \right) + \\ + \,\,{{\alpha }_{p}}\varepsilon {{\varepsilon }_{M}}\frac{{\partial \omega }}{{\partial t}} = {{\varepsilon }_{M}}A(\omega )({{p}_{2}} - {{p}_{1}}), \\ \end{gathered} $(69)
$\begin{gathered} \frac{{\partial {{p}_{2}}}}{{\partial t}} + {{\varepsilon }_{{12}}}\frac{{\partial {{p}_{1}}}}{{\partial t}} + {{b}_{1}}{{k}_{b}}\varepsilon \frac{\partial }{{\partial t}}\left( {\frac{{\partial u}}{{\partial r}} + \frac{u}{r}} \right) - \\ - \,\,\frac{\partial }{{r\partial r}}\left( {{{k}_{2}}({{p}_{2}})r\frac{{\partial {{p}_{2}}}}{{\partial r}}} \right) = A(\omega )({{p}_{1}} - {{p}_{2}}), \\ \end{gathered} $(70)
$\left( {2 + \lambda } \right)\frac{\partial }{{\partial r}}\left( {\frac{1}{r}\frac{{\partial {\mkern 1mu} r{\mkern 1mu} u}}{{\partial r}}} \right) = {{b}_{1}}\frac{{\partial {{p}_{1}}}}{{\partial r}} + {{b}_{1}}{{k}_{b}}\frac{{\partial {{p}_{2}}}}{{\partial r}} + \alpha \frac{{\partial \omega }}{{\partial r}},$(71)
$\frac{{\partial \omega }}{{\partial t}} = \frac{1}{{\tau \beta }}\left\langle {\alpha {{I}_{1}} + {{\alpha }_{J}}J + {{\alpha }_{{p1}}}{{p}_{1}} + {{\alpha }_{{p2}}}{{p}_{2}} - \beta \omega - \gamma } \right\rangle ,$Принимая сделанные выше оценки пороупругих коэффициентов, получим, ${{\varepsilon }_{M}} \approx {{10}^{{ - 2}}}$, $\varepsilon \approx {{10}^{2}}$, ${{\varepsilon }_{{12}}} \approx - 1$, ${{k}_{b}} \approx {{10}^{{ - 1}}}$. В расчетах принималась степенная зависимость обменного члена от параметра поврежденности, $A\left( \omega \right) = A{{\omega }^{2}}$.
РЕЗУЛЬТАТЫ РАСЧЕТОВ
Система дифференциальных уравнений (68)–(71) решалась численно в области 1 < r < 102 на основе полностью неявной схемы. Графики на рис. 4–рис. 6 получены при следующих значениях параметров: ${{\sigma }_{0}} = - 0.0022$,${{p}_{0}} = 0.001$,${{p}_{a}} = 0.0001$, ${{\varepsilon }_{M}} = 0.013$, $\varepsilon = 50$, ${{\varepsilon }_{{12}}} = - 0.19$, ${{b}_{1}} = 0.6$, ${{b}_{2}} = 0.1$, $\lambda = 0.428$, $A = 3 \cdot {{10}^{9}}$, $\alpha = 0.5$, ${{\alpha }_{J}} = 0.5$, ${{\alpha }_{{p1}}} = 0.11$, ${{\alpha }_{{p2}}} = - 0.01$, $\beta = 0.5$, $\gamma = 8 \times {{10}^{{ - 5}}}$. На рис. 4 показана эволюция поврежденности и поровых давлений в матрице и системе трещин. На временах t ~ 1 наблюдается быстрый рост поврежденности в ближайшей (r < 2) окрестности скважины. В дальнейшем при t > 10 поврежденность в этой зоне не испытывает заметных изменений. Под влиянием волны депрессии, распространяющейся от границы скважины r = 1 вглубь пласта, происходит изменение параметров состояния (упругих напряжений и деформаций). Следствием этих факторов является волна разрушения, которая формирует область дренирования матрицы (граница области дренирования приближенно соответствует положению фронта разрушения, ограничивающего поврежденную область, см. рис. 4). Зависимость радиуса границы зоны разрушения d (области дренирования матрицы) от времени показана на рис. 5. Максимальная скорость границы области разрушения наблюдается при малых t. При 102 < t < 103 закон движения границы d(r) аппроксимируется степенной зависимостью. Если принять следующие размерные параметры ${{\mu }_{f}} = {{10}^{3}}$ Па с, $\mu = 4 \times {{10}^{9}}$ Па, $k_{2}^{0} = {{10}^{{ - 12}}}$ м2, получается оценка для характерного времени ${{t}_{0}} \approx 50$ с. Таким образом, согласно рис. 5, радиус границы зоны разрушения достигнет величины $3{{r}_{0}} = 30$ см за 103 с.
Разрушение матрицы приводит к появлению ее проницаемости и перетеканию жидкости в трещины, что приводит к повышению дебита. Факторами, способствующими повышению дебита, являются рост зоны дренирования и увеличение поврежденности. Изменение дебита на единицу длины скважины (${{q(t)} \mathord{\left/ {\vphantom {{q(t)} {{{q}_{r}}}}} \right. \kern-0em} {{{q}_{r}}}}$, где ${{q}_{r}} = {{{{p}_{0}}k_{2}^{2}} \mathord{\left/ {\vphantom {{{{p}_{0}}k_{2}^{2}} {{{\mu }_{f}}}}} \right. \kern-0em} {{{\mu }_{f}}}}$) с течением времени показано на рис. 6 (кривая 1), на больших временах эта зависимость близка к степенной. Для сравнения показан дебит скважины при отсутствии обмена между трещинами и матрицей (кривая 2). В рамках проведенного расчета рост зоны дренирования матрицы не ограничен, однако следует иметь в виду, что максимальная достигаемая поврежденность в волне разрушения быстро падает с увеличением r (см. рис. 4). При достижении параметром поврежденности значений, при которых микротрещины не образуют связанную систему, дальнейший рост зоны дренирования матрицы неизбежно прекратится.
ВЫВОДЫ
В работе представлена модель двойной пористой среды с упругим трещиноватым скелетом с аномально высоким давлением в матрице в изотермическом приближении. Модель согласована с основными принципами механики сплошных сред и термодинамики. Модель включает в себя параметр разрушения, как внутреннюю переменную, характеризующую процесс разрушения в низкопроницаемой матрице вследствие комбинации изменения напряженного состояния и изменения порового давления. Закон образования разрушений сформулирован с учетом неотрицательности диссипации разрушений ${{\delta }_{\omega }}$ для любого механического процесса.
Показано, что аномально высокое давление может привести к развитию разрушений в низкопроницаемой матрице. Подобные условия возникают вследствие более быстрого уменьшения давления в трещинах по сравнению с поровым давлением в матрицах и вследствие изменения полного напряжения в скелете в окрестности скважины. При оценочных значениях параметров модели размер области дренирования первоначально непроницаемой матрицы в результате разрушения в волне депрессии существенно превысил радиус скважины.
Таким образом, продемонстрирована принципиальная возможность вовлечения в разработку низкопроницаемых трещиновато-пористых пластов с АВПД, с использованием только энергии АВПД.
Список литературы
Голф-Рахт Т.Д. Основы нефтепромысловой геологии и разработки трещиноватых коллекторов. М.: Недра. 1986. 608 с.
Извеков О.Я., Кондауров В.И. О рассеянном разрушении пористых материалов с хрупким скелетом // Механика твердого тела. 2010. № 3. С. 182–205.
Извеков О.Я., Кондауров В.И. Модель пористой среды с упругим трещиноватым скелетом // Физика Земли. 2009. № 4. С. 31–42.
Кондауров В.И. Механика и термодинамика насыщенной пористой среды. М.: МФТИ. 2007. 310 с.
Кочарян Г.Г. Геомеханика разломов. М.: ГЕОС. 2016. 424 с.
Мелик-Пашаев В.С., Халимов Э.М., Серегина В.Н. Аномально высокие пластовые давления на нефтяных и газовых месторождениях. М.: Недра. 1983. 181 с.
Мухамедиев Ш.А. Методы локальной реконструкции тектонических напряжений по кинематическим данным: физическая несостоятельность и ложные цели. Ч. I // Физика Земли. 2018. № 6. С. 1–33.
Трусделл К.А. Первоначальный курс рациональной механики сплошных сред. М.: Мир. 1975. 592 с.
Abousleiman Y., Nguyen V. Poromechanics response of inclined wellbore geometry in fractured porous media // J. engineering mechanics. 2005. V. 131. № 11. P. 1170–1183.
Abousleiman Y.N., Tran M.H., Hoang S., Bobko C.P., Ortega A., Ulm F.J. Geomechanics field and laboratory characterization of the Woodford Shale: The next gas play. SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. 2007.
Bai M., Roegiers J.C., Elsworth D. Poromechanical response of fractured-porous rock masses // J. Petroleum Science and Engineering. 1995. V. 13. № (3–4). P. 155–168.
Barati R., Liang J.T. A review of fracturing fluid systems used for hydraulic fracturing of oil and gas wells // J. Applied Polymer Science. 2014. V. 131. № 16.
Barenblatt G.I., Zheltov Y.P., Kochina I.N. 1960. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks // PMM. V. 24. № 5. P. 852–864.
Berryman J.G., Pride S.R. Models for computing geomechanical constants of double-porosity materials from the constituents' properties // J. Geophysical Research: Solid Earth. 2002. V. 107. № B3.
Coussy O. Poromechanics. John Wiley & Sons. 2004.
Craft B.C., Hawkins M.F. Petroleum reservoir engineering. Englewood Cliffs. NJ. 1959. 355 p.
Engelder T., Lacazette A. Natural hydraulic fracturing. Rock joints: Rotterdam, AA Balkema. 1990. P. 35–44.
Eseme E., Urai J.L., Krooss B.M., Littke R. Review of mechanical properties of oil shales: implications for exploitation and basin modeling // Oil Shale. 2007. V. 24. № 2.
Fertl W.H., Chilingarian G.V., Rieke H.H. Abnormal Formation Pressures. American Elsevier Publishing Company. New-York. 1976.
Grytsenko T., Galybin A.N. Numerical analysis of multi-crack large-scale plane problems with adaptive cross approximation and hierarchical matrices //Engineering analysis with boundary elements. 2010. V. 34. № 5. P. 501–510.
Kondaurov V.I., Izvekov O.Y. A model of saturated porous media with elastic brittle skeleton. Proc. Of the 4-the Biot Conference on Poromechanics (Pennsilvania: EStech Publications). 2009. P. 314–20.
Lemaitre J. A Course on Damage Mechanics. Springer-Verlag Berlin Heidelberg. 1996. 228 p.
Li S., George J., Prudy C. Pore-pressure and wellbore-stability prediction to increase drilling efficiency // Journal of Petroleum Technology. 2012. V. 64. № 02. P. 98–101.
Luo X., Vasseur G. Natural hydraulic cracking: numerical model and sensitivity study // Earth and Planetary Science Letters. 2002. V. 201. № 2. P. 431–446.
Magara K. Compaction and fluid migration. Elsevier. 1978.
Mehrabian A., Abousleiman Y.N. Gassmann equations and the constitutive relations for multiple-porosity and multiple-permeability poroelasticity with applications to oil and gas shale. International Journal for Numerical and Analytical Methods in Geomechanics. 2015. V. 39. № 14. P. 1547–1569.
Murakami S. Continuum Damage Mechanics. Springer Netherlands. 2012. 402 p.
Olson J.E. Predicting fracture swarms–The influence of subcritical crack growth and the crack-tip process zone on joint spacing in rock // Geological Society. London. Special Publications. 2004. V. 231. № 1. P. 73–88.
Secor Jr. D.T. Role of Fuid pressure in jointing // American Journal of Science. 1965. V. 263. № 8. P. 633–646.
Snow D.T. Rock fracture spacings, openings, and porosities // J. Soil Mechanics & Foundations Division. 1968. V. 94. № 1. P. 73–92.
Sonnenberg S.A., Pramudito A. Petroleum geology of the giant Elm Coulee field, Williston Basin // AAPG bulletin. 2009. V. 93. № 9. P. 1127–1153.
Thompson J.M., Nobakht M., Anderson D.M. Modeling well performance data from overpressured shale gas reservoirs // Canadian Unconventional Resources and International Petroleum Conference. Society of Petroleum Engineer. 2010. SPE-137755-MS.
Tingay M.R., Morley C.K., Laird A., Limpornpipat O., Krisadasima K., Pabchanda S., Macintyre H.R. Evidence for overpressure generation by kerogen-to-gas maturation in the northern Malay Basin // AAPG bulletin. 2013. V. 97. № 4. P. 639–672.
Warpinski N.R., Mayerhofer M.J., Vincent M.C., Cipolla C.L., Lolon E.P. Stimulating unconventional reservoirs: maximizing network growth while optimizing fracture conductivity // J. Canadian Petroleum Technology. 2009. V. 48. № 10. P. 39–51.
Wu Y.S., Li J., Ding D., Wang C., Di Y. A generalized framework model for the simulation of gas production in unconventional gas reservoirs // SPE Journal. 2014. V. 19. № 05. P. 845–857.
Ye Z., Ghassemi A., Riley S. Fracture properties characterization of shale rocks. In Unconventional Resources Technology Conference. San Antonio. Texas. 1–3 August 2016. Society of Exploration Geophysicists, American Association of Petroleum Geologists, Society of Petroleum Engineers. 2016. 1083–1095 p.
Дополнительные материалы отсутствуют.