Численное моделирование распространения трещины гидроразрыва Текст научной статьи по специальности «Физика»

УДК 539.3

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ РАСПРОСТРАНЕНИЯ ТРЕЩИНЫ

ГИДРОРАЗРЫВА

А. В. Акулич, А. В. Звягин

Одним из эффективных методов повышения производительности нефтедобывающих скважин является искусственное создание в пласте трещин с большой проницаемостью, увеличивающих поверхность фильтрации нефти и как результат — рост ее отдачи в скважине. Впервые задача гидравлического разрыва нефтяных пластов была решена в работах Ю.П. Желтова, С. А. Христиановича [1, 2]. Задача остается актуальной, и ее исследование в разных постановках продолжается. Так, некоторые классы автомодельных решений этой задачи рассмотрены в работах [3-9].

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

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

Постановка задачи. Пусть вязкая несжимаемая жидкость втекает с заданным объемным расходом Qo(t) в центр трещины х = 0 (рис. 1).

1 % 1 У е0«> 1 уф:, /)

. Г" V 4

* 1 и

I

РсЫ

I а0

т

Рис. 1. Геометрическая интерпретация задачи о распространении трещины гидроразрыва

Будем считать текущую ширину трещины ,ш(х,Ь) малой величиной по сравнению с ее текущей полудлиной ¡(Ь), тогда задачу движения жидкости можно рассматривать в одномерной постановке. При этом на выделенный малый элемент жидкости длиной йх действует помимо давления вязкая сила Е йх со стороны стенок трещины. Считая течение достаточно медленным, будем пренебрегать силами инерции жидкости. В этом случае уравнения баланса объема и импульса можно представить в следующем виде:

дги д(тУ) дЬ дх '

дР/

(1)

дх

= - Е,

где Ь — время; V(х, Ь) — скорость жидкости; Р^ — давление в жидкости; Е(х, Ь) — плотность вязких сил, действующих на единицу длины жидкости со стороны стенок. Плотность вязких сил можно принять

такой же, как в задаче о течении вязкой жидкости в подшипнике:

где ц — вязкость жидкости. Мы считаем, что на бесконечности задано давление, т.е. uxx = ayy = —Сто.

Система уравнений (1), (2) незамкнута. Если среда рассматривается как упругая, то в рамках линейной механики трещин и квазистатического приближения уравнения (1) и (2) решаются совместно с уравнениями равновесия упругой среды

даху _ q дсгуХ доуу _ ^

дх ду ' дх ду

и законом Гука:

_ дих_ _ 1 + v г, _ , 1 _ диу _ 1 + V г _ _

£хх — Q — Lv1 1У)<тхх иауу\> £УУ — Q — Lv1 1У)<ТУУ U<Txx\ ,

У (4)

1 ( дих диу \ 1 + v

£ху =£ух = ^ [-¿^ + — <тху.

Здесь x, y — координаты; axx, axy, ayy — компоненты тензора напряжений; exx, exy, eyy — компоненты тензора деформаций; ux, uy — компоненты вектора перемещений; E — модуль Юнга; v — коэффициент Пуассона.

Уравнения (1), (2) и (3), (4) связаны граничными условиями на берегах трещины (y = 0±, \x\ < l). Введем следующие обозначения для параметров на берегах трещины:

а± = ai3(x, 0±), u± = Ui(x, 0±), w(x,t) = u+(x, 0,t) — u-(x, 0,t).

Тогда условия на берегах принимают форму

y = 0±, \x\ <l(t), a±y = 0, a±y = —Pf(x,t). (5)

Здесь необходимо заметить, что точным граничным условием было бы условие a±y = —F, а не a±y = 0, но мы вслед за [3-6, 8, 9] считаем, что сила F по модулю много меньше давления Pf (x,t), и пренебрегаем ею.

Кроме того, необходимо добавить условия втекания жидкости и критерий разрушения среды. Пусть Vo(t) — скорость жидкости, wo(t) — раскрытие трещины в точке с координатой x = 0, а Q = Qo/2 —

половина расхода жидкости. Тогда в силу симметрии задачи

Vb(i) = (6)

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

К1 = 7ш1(Р/(е'" ^ \[Щ= Klc~ (7)

-1 *

Здесь К — коэффициент интенсивности напряжений, К1С — его критическая величина.

В рассматриваемой постановке задачи положение переднего фронта жидкости считается неизвестным и определяется координатой х = Ь(Ь) (рис. 1). Из закона сохранения объема для закачанной жидкости получаем соотношение

ь г

/ w(x,t) (1х = Q(t) (И, (8)

где Ь(Ь) — положение фронта жидкости в трещине. Давление жидкости на этом фронте должно быть равно нулю:

х = щ, Рг (ь,г) = о. (9)

Кроме того, на бесконечности в упругой среде должны быть выполнены следующие граничные условия:

х2 + У2 ^ Ж, аXX = °уу = -&0, &ху = 0. (10)

Таким образом, задача сводится к решению уравнений (1)-(4) с учетом условий (5)-(10).

Метод построения решения. Анализ задачи показывает, что в начальный момент времени, когда длина трещины равна нулю, давление в трещине должно быть бесконечным, поэтому при численном решении необходимо выбрать начальное приближение в некоторый момент времени Ь = ¿о, отличный от нуля, чтобы начать процесс вычислений. В качестве таких начальных условий было взято точное решение, полученное в [8]:

Q3/4 У0 =-*

^^(f) 1/4 ' «*(*) = т7/Ч1о - х)^ (!§-)1/3, Мх)=р> -

6^2 /(£')У^о\1/3 г, Е „ , , , ,

где р = do H----- ; Ь--Vo, WqIx) и ро{х) — скорость жидкости, профиль раскрытия

п \ 3 l0 J 1 — v2

трещины и профиль давления в начальный момент времени соответственно. Расчет начинался с длины lo = 1 м, что соответствует моменту времени to = 0,63 с.

Задача определения перемещений и напряжений для упругой среды в данной работе решается методом разрывных смещений [10, 11], который является очень эффективным для решения задач механики трещин. Уравнение сохранения массы (1) линейно относительно x и может быть легко проинтегрировано. Тогда для скорости V(x,t) в каждый момент времени получим выражение

x

c — f W(x, t) dx

Принимая в расчет граничное условие в зоне закачки жидкости (6), будем иметь

Q

x = 0, V (0,t) =

wo(t) '

откуда определяется с = Q — константа интегрирования в (11). Для расчета скорости жидкости внутри трещины предположим, что раскрытие трещины представлено в следующей форме [7]:

w(x,t)=w0(t)f^щj =и>о№, (12)

х

где £ = ——- — безразмерная координата. Интеграл в выражении для скорости (11) примет вид ¡(Ь)

X ? ? ?

У й(х, Ь) йх = ^¡1/ (С) йС - С/ '(С) йС = №о I + wоi)J С/ (С) йС - wоС/(01. (13)

0 0 0 о

Из закона сохранения объема следует, что для постоянного расхода Q

1

/ Ж) ¿е = оь, гЬ01 + и}01 = у-^-• (14)

о / / (С) йС

0

Подставляя (14) в (13) и затем в (12), мы можем найти распределение скорости как функцию расхода жидкости Q и скорости конца трещины I, а также зависимость профиля скорости жидкости от формы раскрытия трещины:

(

V (£,t) =

Q

w

1 -

/ f (t) d^

0

V

I f (t) d£

0

+1

(15)

Численный метод. Решая задачу, мы будем использовать факт, что задача теории упругости эллиптического типа, и поэтому ее решение непрерывно зависит от параметров. В данной постановке время играет роль параметра, поскольку все уравнения, за исключением уравнения неразрывности, приняты статическими. При этом следует иметь в виду, что система уравнений движения жидкости (1), (2) и система уравнений статической теории упругости (3), (4) должны решаться совместно для каждого момента времени. В наших вычислениях мы применяем метод последовательных приближений. В качестве начального приближения в каждый момент времени используется автомодельное решение [7]. С целью увеличения точности расчетов при сохранении общего количества граничных элементов мы прибегаем к неравномерному разбиению, самые малые граничные элементы будут у конца трещины. Для этого мы при замене £ = sin р используем переменный шаг по р, где 0 < р < п/2.

Для выбранного разбиения по длине трещины в безразмерных координатах получаем следующую запись — дискретную форму уравнения (15):

T/fc tkf , Q (л Ik((N)

v =tL+^(l-im

(16)

N

к

где Ik (£N) = Y,

wk l+wk d£k k k _ ^

m=1 m=1

w

k-1

+ wk d£k

— , к w0

номер узла. Чтобы найти wk, мы ис-

2 адо 1 I2

т=1

пользуем в качестве первого приближения определенную на предыдущем временном шаге величину ,шк. Для интегрирования уравнения закона сохранения импульса (1) оно переписывается в следующем дис-

кретном виде:

pk = pk-1 - 12v

Vk-1 + Vk 2

(17)

Задача теории упругости (3), (4) решается методом разрывных смещений.

Метод разрывных смещений основан на представлении решения задачи теории упругости в виде суммы базисных аналитических решений с неизвестными коэффициентами. Сами коэффициенты определяются из условий выполнения граничных условий на дискретном множестве точек границы. В качестве базисных решений берутся решения задач о постоянном разрыве перемещения на конечном отрезке в упругой плоскости (граничном элементе). Трещина разбивается на N отрезков, и каждый из них используется в качестве граничного элемента. Зная аналитическое решение для единичного разрыва перемещений, находим численное решение задачи, суммируя эффекты всех N элементов. Подробное описание метода можно найти в [7].

Общая схема решения задачи для каждого шага по времени сводится к последовательным итерациям выбранного начального приближения с целью минимизации погрешности. Итерационными параметрами задачи для фиксированного значения времени Ь являются скорость ¡(Ь) трещины и скорость Ь(Ь) переднего фронта жидкости. При известном решении на предыдущем шаге по времени, по которому строится начальное приближение, и заданных скоростях ¡(Ь), Ь(Ь) уравнения (16), (17), а также уравнения (3), (4) с граничными условиями (5)—(10) интегрируются. Два вышеназванных параметра находятся из условия выполнения критерия разрушения в конце трещины и условия совпадения объема жидкости в трещине с закачанным объемом:

\Ki - KIc\ <£1,

L(t) t

j w(x, t) d£ - j Q(t) dt 00

< £2,

где е\, £2 — соответствующие погрешности. Данная схема решения была реализована в виде программы. С целью проверки ее работы проведены серия расчетов и их сравнение с результатами других авторов.

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

1) трещина полностью заполнена закачиваемой жидкостью (иными словами, фронт жидкости совпадает с концом трещины);

2) упругая среда не имеет прочности на разрушение (критический коэффициент интенсивности напряжений равен нулю).

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

При решении автомодельной задачи вводились следующие величины:

Т = —^ — временной масштаб, Ь' = \ZQoT — масштаб длины. (18)

Е

Безразмерные координата С и время т определялись следующим образом:

{ = т = ± (19)

Аналогично вводились безразмерная длина трещины х(т), давление П(С,т) и раскрытие трещины 0(С,Ь):

х = V' п = » = £ <2»>

Автомодельное решение для раскрытия и давления имеет форму

п(С,т ) = Хт 1/3Щ (С), П(С, т ) = т-1/3 п (С), (21)

где функция П (С) представляется суммой двойного ряда и имеет достаточно сложный вид [4]. Согласно [4], Х ^ 0,629.

Для проведения сравнения расход жидкости, параметры упругой среды и жидкости выбирались такими, как и в [4], и строились графики различных величин в зависимости от времени или от координаты при фиксированном времени. Результаты численных расчетов приводились к форме (18)-(22) и сравнивались с аналитическими.

В численной модели использовались следующие входные данные. Параметры упругой среды: Е = 25 ГПа; V = 0,15; а0 = 50 МПа; Q0 =4 • 10-3 м2/с. Вязкость жидкости ц = 8,53 • 10-7 МПа • с. Расчеты проводились для следующих значений прочности среды: К1С = 0,001; 0,1; 1; 5; 10 (значения взяты в

единицах МПа • м1/2). В безразмерном виде прочность определялась по формуле

*=(22)

То есть в вычислениях использовались следующие величины безразмерной прочности: К = 3,5 • 10-5; 3,5 • 10-3; 3,5 • 10-2; 1,7 • 10-1 и 3,5 • 10-1.

На рис. 2 приведены графики давления в центре трещины и раскрытия трещины в ее центре в зависимости от времени, а также зависимость длины трещины от времени. Сравнение численного и автомодельного решений для величин К вплоть до 1,7 • 10-1 указывает на очень хорошее согласование этих решений.

На рис. 3, а графики П^(С) и (С) приведены для разных фиксированных моментов времени при значении К = 3,5 • 10-5. Существенное различие между аналитическим и численным решениями наблюдается для малых моментов времени, когда прочность среды и наличие отставания жидкости от конца трещины оказывают заметное влияние на процесс. Интересно, что при большой длине трещины автомодельное и численное решения снова совпадают. Для очень малой прочности решения фактически не отличаются, что является хорошей тестовой проверкой для используемой оригинальной программы. Результаты для К = 3,5 • 10-3; 3,5 • 10-2 и 1,7 • 10-1 не содержат существенных отличий от случая К = 3,5 • 10-5.

(И

1 х К=?,5- КГ5 III

+ 3,5 -10"2 □ 1.7-10"'

ж 3,5 чо-!

И

у

-В- 1 I

к-— —ft—I— *—*—в— 4—й--я

0 0,2 0,4 0,6 0,8 1 1,2

Рис. 2. Решение в точке закачки жидкости. Зависимость от безразмерного времени т =

т ■ Ю-12 = — 10~12 безразмерного давления П = П • 104 = ^ —- 104 (кривая /); без-T E

размерной длины трещины \ = X ' Ю-7 = — Ю-7 (кривая II); безразмерного раскры-

L

тия Q = Q • 10~3 = — 10~3(кривая III). Сплошная линия — аналитическое решение L

а

:! i с ) X 1 X т=1,9- 10ч д 1J-1010

/ X к + ж 0,! О12

ж"—-- X

Л X

II \

а- Г-- X Л V

X 1 х T=l,9-104 й 1,3 *1В10

I + ж О1? „12

Щ,

II

* X | д V \ + V

-^-iSst

0 0,2 0,4 0,6 0,К 1 0 0,2 0А 0,6 0,К 1

Рис. 3. Решения для случая безразмерной прочности К = 3,5 • 10~5 (а) и К = 3,5 • 10-1 (б). Зависи-

X "

мость от безразмерной координаты £ = — усеченного безразмерного давления П^(^) = П^(^) • 10 =

г1/3П(^, г) • 10 (кривая I); усеченного безразмерного раскрытия трещины = —

Л

(кривая II). Сплошная линия — аналитическое решение

Соответствующие графики для большого значения прочности К = 3,5 • 10_1 показаны на рис. 3, б. В этом случае численные результаты обнаруживают значительное расхождение с полуаналитическим решением. Это объясняется тем, что здесь аналитическое решение, полученное при отсутствии прочности, неприменимо.

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

СПИСОК ЛИТЕРАТУРЫ

1. Khristianovich S.A., Zeltov Y.P. Formation of vertical fractures by means of highly viscous liquid // Proc. Forth World Pet. Congress, Rome. Vol. 2. 1995. 579-586.

2. Христианович С.А., Желтое Ю.П. О гидравлическом разрыве нефтеносного пласта // Изв. АН СССР. ОТН. 1955. № 5. 3-41.

3. Perkins N.K., Kern L.R. Widths of hydraulic fractures // Paper SPE 89. J. Petroleum Technol. 1961. 13, N 9. 937-949.

4. Adachi J. I., Detournay E. Self-similar solution of a plane-strain fracture driven by a power-law fluid // Int. J. Numer. Anal. Methods Geomech. 2002. 26. 579-604.

5. Spence D.A., Sharp P. Self-similar solution for elastohydrodynamic cavity flow // Proc. Roy. Soc. London. Ser. A. 1985. 400. 289-313.

6. Ивашнев О.Е., Смирнов Н.Н. Формирование трещины гидроразрыва в пористой среде // Вестн. Моск. ун-та. Матем. Механ. 2003. № 6. 28-36.

7. Звягин А.В. Движение вязкой жидкости в канале с упругими стенками // Вестн. Моск. ун-та. Матем. Механ. 2005. № 1. 50-54.

8. Desroches J., Detournay E, Lenoach B, Papanastasiou P., Pearson J.R.A., Thiercelin M., Cheng A.H.-D. The crack tip region in hydraulic fracturing // Proc. Roy. Soc. London. Ser. A. 1994. 447. 39-48.

9. Carbonell R., Desroches J., Detournay E. A comparison between a semi-analytical and a numerical solution of a two-dimensional hydraulic fracture // Int. J. Solids Structures. 1999. 36. 4869-4888.

10. Crouch S.L., Starfield A.M. Boundary element methods in solid mechanics. London; Boston: Allen & Unwin, 1983.

11. Богданов А.В., Звягин А.В., Тьерсилен М. Взаимное влияние системы трещин на коэффициент интенсивности напряжений // Вестн. Моск. ун-та. Матем. Механ. 2004. № 6. 44-49

Поступила в редакцию 10.07.2006

УДК 534.222

ВЛИЯНИЕ ДИФФУЗИОННЫХ И КИНЕТИЧЕСКИХ ХАРАКТЕРИСТИК НА ПРОЦЕСС ГОРЕНИЯ ЧАСТИЦЫ УГЛЕРОДА

Н. Н. Смирнов, В. Н. Янушкевич

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

1. Решение задачи горения частицы углерода с учетом выхода летучих компонентов. Рассмотрим следующие реакции, протекающие в конденсированной фазе:

Ь(р) ^ Ь,

2С + 02 ^ 2СО, (1)

С + СО2 ^ 2С0. (2)

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

Настоящая модель включает основополагающие части моделей Нуссельта и Бурке-Шумана [1, 2].

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

Окисление углерода (1), (2) определяется химической кинетикой реакций на поверхности и диффузией реагентов к поверхности. Выделение летучих веществ разбавляет смесь, влияя таким образом на скорость реакций на поверхности.