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

© И.В. Моренко, Ю.А. Волков УДК 533.2

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

И.В. Моренко1, Ю.А. Волков2

Обособленное структурное подразделение ФИЦ Институт механики и машиностроения КазНЦ РАН, г. Казань, Россия 2ЦСМРнефть, г. Казань, Россия

то гепко @1тт. knc.ru

Резюме: Рассмотрены гидродинамические процессы в обсадной колонне вертикальной скважины при имплозии. Математическое моделирование движения двухфазной среды со свободной поверхностью основано на методе объема жидкости. Для замыкания осредненных уравнений Рейнольдса используется зональная модель Ментера SST. Численные расчеты выполнены методом конечных объемов с помощью открытой интегрированной платформы OpenFOAM. Установлено, что при заданных входных параметрах период депрессии составляет 1 с, а максимальное давление в контрольной точке превосходит гидростатическое в 1,8 раза.

Ключевые слова: численное моделирование, двухфазная среда, метод объема жидкости, имплозия.

Вв1: 10.30724/1998-9903-2018-20- 1-2 -128-137

CALCULATION OF NON-ISOTHERMAL MOTION OF A TWO-PHASE COMPRESSIBLE MEDIUM AT IMPLOUSION

I.V. Morenko1, Y.A. Volkov2

1 IME KazSC RAS, Kazan, Russia 2 Oil field development methods improvement center ltd, Kazan, Russia

morenko @imm. knc. ru

Abstract: Hydrodynamic processes in a vertical well casing during implosion are considered. Mathematical modeling of the motion of a two-phase medium with a free surface is based on the volume of fluid method. Unsteady Reynolds-Averaged Navier-Stokes equations are supplemented by the Menter's Shear Stress Transport turbulence model. Numerical calculations were performed by the finite volume method using the open integrated platform OpenFOAM. It was found that the depression period is 1 s, and the maximum pressure at the control point is 1.8 times higher than the hydrostatic pressure at the given input parameters.

Keywords: Numerical simulation, two-phase medium, volume of fluid method, implosion. Введение

Течение жидкости со свободной поверхностью имеет место в природе. Например, это: течение рек, набегание волн на берег, дамбу и опоры моста [1]; плескание жидкости в технических устройствах, в топливных баках самолетов и космических ракет, в цистернах

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

В нефтедобывающей отрасли имплозионная технология широко используется для улучшения фильтрационных характеристик призабойной зоны пласта. С этой целью герметически закрытая заполненная воздухом имплозионная камера низкого давления (как правило - цилиндрической формы) опускается в нижнюю часть (на уровень пласта) заполненной жидкостью скважины. Снизу она закрыта мембраной. Когда мембрана «мгновенно» устраняется, камера заполняется скважинной жидкостью. При этом воздух, находившийся в камере, сильно сжимается. Имплозионный эффект на первой стадии воздействия заключается в мгновенном снижении давления в призабойной зоне пласта, что должно способствовать выносу загрязнений из него в случае наличия перфорации в обсадной колонне. Вторая стадия имплозионного воздействия - это гидравлический удар падающего столба скважинной жидкости, способный создать в нижней части скважины и, в частности, в призабойной зоне пласта давление, значительно превышающее пластовое и способствующее, при наличии перфорации, образованию в пласте трещин [2; 3].

До сих пор для расчета имплозионного процесса в стволе скважины, вскрывающей нефтяной пласт, использовались полуэмпирические математические модели и аналитические решения [2-5]. Однако развитие численных методов и наличие высокопроизводительных компьютеров позволяет рассчитывать процессы, протекающие при имплозии с использованием математических моделей, основанных на законах сохранения, и решать дифференциальные уравнения в частных производных.

К настоящему времени для математического моделирования движения жидкости со свободной поверхностью в пакетах программ с открытым исходным кодом реализованы следующие методы [6; 7]: Volume of Fluid (метод объема жидкости), Smoothed Particle Hydrodynamics (метод сглаженных частиц), Particle Finite Element Method (метод конечных элементов с частицами).

В данной работе для расчета нестационарного неизотермического движения жидкости со свободной поверхностью используется метод объема жидкости (Volume of Fluid) [8], согласно которому течение среды моделируется единым набором уравнений движения, энергии для всех фаз.

Постановка задачи

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

Расчетная область Q представляет собой цилиндрический сектор с угловым размером 0 = 5 ° (рис. 1). Центр декартовой системы координат XjX2X3 разместим так,

чтобы ось Ox3 совпадала с осью обсадной колонны с высотой hg , внутренним радиусом rt. Газом заполнена подобласть Qg, а также имплозионная камера (область низкого давления) а, в подобласти а находится жидкость. Границу обсадной колонны скважины обозначим через Гs, имплозионной камеры - Г^, мембрану имплозионной камеры - Гm . Длина имплозионной камеры hk2 обычно сопоставима с толщиной пласта.

Рассматривается двухфазная среда, которая состоит из жидкой и газовой фаз. Обозначим через а объемную долю жидкой фазы. Если контрольный объем занят жидкостью, то а = 1, если газом - а = 0, если через контрольный объем проходит фазовая граница, то 0 < а < 1.

Рис. 1. Схема расчетной области

Плотность р , динамический коэффициент вязкости ц, удельная теплоёмкость при постоянном давлении ер, коэффициент теплопроводности X среды рассчитываются по соответствующим значениям жидкости (индекс I) и газа (индекс g) с использованием маркерной функции а:

р=р1а+рё (1-а) , ц=ц1а+цё (1-а),

S =(^ )г а + (ср ) (1 - а) , к=

hХ g

ha+hg (1-а)

(1)

Запишем уравнения состояния газа и жидкости

рg = pRgT ; Pi = Pl0 + PliRlT

здесь Rg - универсальная газовая постоянная; R = const ; р/0 = const.

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

до d(pU )

- + -

dt дх.

■ = 0.

Уравнение количества движения д(рщ ) d(pu1uj )

dt

- + -

dx.

= Pg

dp d dx, dx,

(

Veff

\

dUi du^

dx, dx, V J ' j

- 3 «J

fdu± \dxk j

+oxVa.

(2)

(3)

"1 ""I 1

Здесь и - компонента осредненного вектора скорости и; / - время; р - давление; = ц + - эффективная вязкость; цт - турбулентная, или вихревая вязкость; 5у -

символ Кронекера; g^ - компоненты ускорения свободного падения; /с = -У- п - кривизна

свободной поверхности; п =Уа/|Уа| - единичная нормаль к поверхности раздела фаз; а -

коэффициент поверхностного натяжения.

Для определения положения свободной поверхности записывается уравнение конвективного переноса

^ + У-(аы ^-((1 - а) аиг ) = 0 . (4)

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

иг = и1 -и8 .

Уравнение сохранения энергии

д(pcpT) д

дt

-(рсрТы1)--

дх V^ 1 Р^

дх.

дТ дх,

= 0.

к У

(5)

ср Дт

здесь Т - температура; 1 ~ = 1 + ——; Ргг = 0.85 - турбулентное число Прандтля.

Рг,

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

Моделирование турбулентности осуществляется на основе хорошо апробированной [9] зональной модели Ментера [10]. Модель Ментера SST включает уравнения для транспорта кинетической энергии турбулентности к и скорости диссипации турбулентной энергии ш:

д(рк ) д(ры}к )

д/

дх.

= Р - в ршк +

_д_

дх

(ц+акДт)

дк дх

(6)

д(рш) , дК ») = Тга. , - рюрй2 + А

д/

дх.

дх.

/ ч дш

+ 2 (1-Л ) ра £ £. (7)

ш дх^ дх^

х1 Vт "х1

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

Р = шт

т (Р, 20в *шк), Р =

диг

дх.

тензор вязких напряжений т^ = 2дТБу - — рк5у

тензор скоростей деформации =

Г ды, ды!Л —- + —-

5х, дх, V 1 '

. Коэффициент турбулентной вязкости

Дт =

ра1к

шах [о^ш, ОЕ2 ]

щ = -1 2

( ды, ды- )

дх, дх, V - г у

^ = 1апЬ (а^ ),

(

aгg2 = шах

4к 500у

в шй ^ ш

Функция переключения ^ имеет вид ^ = 1апЬ) .

aгg1 = шт

шах

•ч/к 500V

в шй й ш

4аш2 к СОкшй:

(

СГ^ = шах

Л

2аш21 — до-20 ш дх,- дх,

V 11 У

й -

расстояние от узла сетки до стенки.

г1

Константы уш, Рш , рассчитываются по формулам: Рю = +(1 -F1)Рю2 ;

о *2 о *2

к -ГТгт ./"i . _ в1 °ю1К „ _Р2 °ю2К

°га = F1°ra1 + 0 - F1) °и2 ; Yl = -¡T--, Y2

Р* Тв7 2 р* # ■

Константы модели он = 0,85 ; 2 = 1,0 ; ош1 = 0,5 ; ош2 = 0,856 ; Р1 = 0,075 ; Р2 = 0,0828; р* = 0,09 ; к* = 0,41; а = 0,31.

Граничные и начальные условия задаются следующим образом. На поверхностях Г = Гж и Г^ обсадной колонны скважины Гж и имплозионной камеры Гк записывается

I А Ф

условие прилипания, температура и = 0, —

1 дп

= 0, ^ дп

= 0, Г|г = Tw, ¿г = 0 ,

г

, 6v

ю = 10-- , здесь Р1 = 0,075 , d - пристеночный шаг, п - единичная нормаль.

1 Ml2

На боковых границах 1side ставятся условия симметрии.

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

Решение поставленной задачи (1)-(7) осуществляется на гексаэдральной сетке методом конечных объемов с помощью открытой интегрируемой платформы для численного моделирования задач механики сплошных сред OpenFOAM [11] с лицензией GNU GPL. Общее число узлов расчетной сетки составляет порядка 105. Сгущение сетки выполнено в областях больших градиентов параметров течения, в частности, в области камеры и вблизи межфазных границ. Дискретные значения скорости и давления рассчитываются в центрах ячеек сетки. Для решения системы линейных уравнений для давления используется алгебраический многосеточный метод, для скорости - метод бисопряженных градиентов. Шаг по времени выбирается таким образом, чтобы число Куранта не превосходило 0,03.

Теплофизические свойства газа и жидкости

Задаются следующие теплофизические свойства газа: коэффициент кинематической вязкости vg = 1,46-105 м2/с, удельная теплоёмкость при постоянном давлении cpg =1006 Дж/(кгК), коэффициент теплопроводности Xg = 0.024 Вт/(м К), а также свойства жидкости: плотность р/0 = 103 кг/м3; коэффициент кинематической вязкости V = 10-6 м2/с; удельная теплоёмкость при постоянном давлении c ¡ =4182 Дж/(кгК); коэффициент теплопроводности X¡ = 0,6 Вт/(мК). Кроме того, задаются Rg = 8,314 Дж/(мольК) - универсальная газовая постоянная, R = 3000 Дж/(мольК), о = 0,07 Н/м - коэффициент поверхностного натяжения.

Тестирование программного комплекса

С целью верификации программного комплекса решается классическая тестовая задача обрушения столба жидкости. Рассмотрим прямоугольный резервуар высотой 0,584 м, длиной 0,584 м, в котором находится столб жидкости высотой b =0,292 м, длиной а =0,146 м в гидростатическом равновесии. Столб жидкости ограничен вертикальной стенкой. В начальный момент времени стенка удаляется, а столб жидкости под действием силы тяжести начинает обрушаться, формируется волна (рис. 2). В ходе численного эксперимента определяется положение фронта волны с течением времени. Для этого обозначим через х* = х/а безразмерное положение фронта волны, t* = t^ng/а -

г

безразмерное время, п = Ъ/а. Сопоставление расчетных значений положения фронта волны и экспериментальных данных [12] показывает корректность результатов, полученных с использованием метода УОЕ (рис. 3).

Рис. 2. Расчетное положение свободной поверхности в моменты времени а) ^ = 0 ; Ь) /* = 1,16 ; с) /* = 2,32; ф) /* = 3,48

Рис. 3. Положение фронта воды в зависимости от времени: __- численный расчет, • - эксперимент [12]

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

Численное моделирование имплозионного процесса выполнено для случая, когда обсадная колонна скважины имеет внутренний радиус г = 0,12 м, высоту = 2000 м.

Высота столба скважинной жидкости 1800 м. Длина имплозионной камеры составляет 50 м, внутренний радиус камеры г1 = 0,065 м, внешний радиус г2 = 0,07 м. Величина зумпфа 2,5 м (часть скважины под дном камеры). Обсадная колонна и пласт - не перфорированы.

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

Распределение давления в радиальном сечении обсадной колонны скважины в первые 2,4 с иллюстрирует рис 4. В начальный момент времени мембрана имплозионной камеры закрыта и система находится в гидростатическом равновесии. С открытием мембраны начинается первый этап имплозионного воздействия - поглощение камерой воды и создание депрессии в скважине. Сначала давление в скважине падает лишь на уровне мембраны. Затем область пониженного давления распространяется от уровня мембраны вниз (в зумпф) и вверх по колонне. Одновременно наблюдается падение давления в направлении от устья (межфазная граница «жидкость -газ») ко дну скважины. В момент времени приблизительно 0,95 с давление по всей длине скважины становится низким. После этого начинается второй этап первого цикла имплозионного воздействия -падение столба жидкости. Волна высокого давления, которое значительно превышает пластовое, начинает распространяться со дна имплозионной камеры вниз (рис. 4 t = 1,0 с), затем вверх по колонне.

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

0.01с 0.2с 0.4с 0.6с 0.8с 1.0с 1.2с 1.4с 1.6с 1.8с 2.0с 2.2с 2.4с

rJ

Рис. 4. Поля давления в радиальном сечении обсадной колонны скважины в разные моменты времени. Для наглядности иллюстрации приведены с разным масштабом по

осям абсцисс и ординат

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

В начальный момент времени давление в точке Р составляет р0 = 1,8 • 107 Па. После открытия мембраны величина давления резко падает в течение 0,004 с. Период депрессии длится около 1 с. При этом давление в точке Р в среднем составляет от 3,3 -106 до 9,4 -106 Па. На данном этапе имплозионного процесса низкое давление р < р0 будет способствовать очистке призабойной зоны.

Рис. 5. Изменение давления (Па) в контрольной точке Р с течением времени (сплошная линия); р0 (пунктирная линия)

Скачок давления происходит в момент времени t«1 с. Рост давления происходит за счет падающего столба жидкости. Начинается вторая стадия первого цикла имплозионного воздействия, продолжительностью около 2 с. Следует отметить, что момент времени, когда жидкость заполняет камеру, соответствует максимальному

7

давлению р = 3,27-10' Па, которое достигается в контрольной точке. Максимальное давление в точке Р при этом превышает гидростатическое давление в 1,8 раза. Кроме того, гидродинамическая сила, действующая на камеру, достигает своего максимума Е = 5,76 • 108 Па. Вторая стадия имплозионного воздействия, характеризующаяся скачком давления, на практике будет сопровождаться образованием трещин в нефтяном пласте.

Далее период высокого давления (р > р0) сменяется низким (р < р0). Начинается второй цикл имплозионного воздействия. Численный расчет показал, что колебания уровня жидкости, поля давления в расчетной области, а также силы, действующей на имплозионную камеру, затухают с течением времени пока система не придет в гидростатическое равновесие.

Заключение

Метод объема жидкости, реализованный на интегрированной платформе ОргпЕОАМ, является эффективным инструментом для моделирования динамики газожидкостной смеси со свободной границей. Показано, что математическое и численное моделирование процесса движения среды после открытия мембраны имплозионной камеры, находящейся в заполненной водой скважине, позволяет с помощью функционала ОргпЕОАМ, исследовать основные этапы имплозионного воздействия, определять их продолжительность, а также в деталях анализировать механизмы процессов, протекающих в системе «скважина-водяной столб со свободной поверхностью на устье-имплозионная камера на уровне пласта».

Литература

1. Vassilevski Yu.V., Nikitin K.D., Olshanskii M.A., Terekhov K.M. CFD technology for 3D simulation of large-scale hydrodynamic events and disasters // Russ. J. Numer. Anal. Math. Modelling, 2012. Vol. 27, No. 4. P. 399-412.

2. Волков Ю.А., Конюхов В.М., Костерин А.В., Чекалин А.Н. Математическое моделирование имплозионного воздействия на пласт. Казань: «Плутон», 2004. 78 с.

3. Попов А.А. Имплозия в процессах нефтедобычи. М.: Недра, 1996. 192 с.

4. Косолапов А.Ф., Корженевский А.Г., Корженевская Т.А. Энергетические, гидродинамические и волновые характеристики имплозионного источника ДСК-1 // Научно технический вестник «Каротажник». 2015. № 247. С. 92-100.

5. Конюхов В.М., Краснов С.В., Конюхов И.В. Компьютерное моделирование имплозионных процессов в системе «нефтяной пласт - скважина - подвижная имплозионная камера» // Научно-технический журнал «Автоматизация, телемеханизация и связь в нефтяной промышленности». 2016. № 5. С. 32-37.

6. Давыдова Е.В., Корчагова В.Н. Свободное программное обеспечение для моделирования жидкости со свободной поверхностью // Труды института системного программирования РАН. 2016. Т. 28. № 1. С. 243-258.

7. Шаланин В.А. Лагранжевы методы моделирования потоков со свободной поверхностью // Молодой ученый. 2016. № 2. С. 261-263.

8. Hirt C.W., Nichols B.D. Volume of fluid (VOF) method for the dynamics of free boundaries // Journal of Computational Physics. 1981. Vol. 39, N. 1. P. 201-225.

9. Моренко И.В., Федяев В.Л. Влияние турбулентности потока вязкой жидкости на гидродинамические характеристики и теплообмен обтекаемых тел // Известия вузов. Проблемы энергетики. 2010. № 7-8. С. 36-45.

10. Мenter F.R. Two-equation eddy viscosity turbulence models for engineering applications // AIAA J. 1994. V. 32. № 8. P. 1598-1605.

11. OpenFOAM - The open source CFD toolbox [Электронный ресурс]. URL: http://www.openfoam.com.

12. Martin J.C., Moyce W.J. An experimental study of the collapse of liquid columns on a rigid horizontal plane // Phil. Trans. Roy. Soc. London (Philosophical transactions of the royal society a mathematical, physical and engineering sciences). 1952. Vol. 244. No. 882. Pp. 312-324. DOI: 10.1098/rsta.1952.0006.

Авторы публикации

Моренко Ирина Вениаминовна - канд. техн. наук, старший научный сотрудник Института механики и машиностроения, обособленное структурное подразделение ФИЦ Институт механики и машиностроения КазНЦ РАН, г. Казань, Россия

Волков Юрий Андреевич - канд. физ.-мат. наук, генеральный директор OOO «ЦСМРнефть».

References

1. Vassilevski Yu.V., Nikitin K.D., Olshanskii M.A., Terekhov K.M. CFD technology for 3D simulation of large-scale hydrodynamic events and disasters // Russ. J. Numer. Anal. Math. Modelling, 2012. Vol. 27, No. 4. P. 399-412.

2. Volkov Ju.A., Konjuhov V.M., Kosterin A.V., Chekalin A.N. Matematicheskoe modelirovanie implozionnogo vozdejstvija na plast. Kazan': «Pluton», 2004. 78 s.

3. Popov A.A. Implozija v processah neftedobychi. M.: Nedra, 1996. 192 c.

4. Kosolapov A.F., Korzhenevskij A.G., Korzhenevskaja T.A. Jenergeticheskie, gidrodinamicheskie i volnovye harakteristiki implozionnogo istochnika DSK-1 // Nauchno tehnicheskij

npoöneMbi энергетики, 2018, moM 20, № 1—2

vestnik «Karotazhnik». 2015. № 247. P. 92-100.

5. Konjuhov V.M., Krasnov S.V., Konjuhov I.V. Komp'jutemoe modelirovanie implozionnyh processov v sisteme «neftjanoj plast - skvazhina - podvizhnaja implozionnaja kamera» // Nauchno-tehnicheskij zhurnal «Avtomatizacija, telemehanizacija i svjaz' v neftjanoj promyshlennosti». 2016. № 5. P. 3237.

6. Davydova E.V., Korchagova V.N. Svobodnoe programmnoe obespechenie dlja modelirovanija zhidkosti so svobodnoj poverhnost'ju // Trudy instituta sistemnogo programmirovanija RAN. 2016. T. 28. № 1. P. 243—258.

7. Shalanin V.A. Lagranzhevy metody modelirovanija potokov so svobodnoj poverhnost'ju // Molodoj uchenyj. 2016. № 2. P. 261—263.

8. Hirt C.W., Nichols B.D. Volume of fluid (VOF) method for the dynamics of free boundaries // Journal of Computational Physics. 1981. Vol. 39, N. 1. P. 201-225.

9. Morenko I.V., Fedjaev V.L. Vlijanie turbulentnosti potoka vjazkoj zhidkosti na gidrodinamicheskie harakteristiki i teploobmen obtekaemyh tel // Izvestija vuzov. Problemy jenergetiki. 2010. № 7-8. P. 36-45.

10. Menter F.R. Two-equation eddy viscosity turbulence models for engineering applications // AIAA J. 1994. V. 32. № 8. P. 1598-1605.

11. OpenFOAM - The open source CFD toolbox [Jelektronnyj resurs]. URL: http://www. openfoam. com.

12. Martin J.C., Moyce W.J. An experimental study of the collapse of liquid columns on a rigid horizontal plane // Phil. Trans. Roy. Soc. London (Philosophical transactions of the royal society a mathematical, physical and engineering sciences). 1952. Vol. 244. No. 882. Pp. 312-324. DOI: 10.1098/rsta.1952.0006.

Authors of the publication

Irina V Morenko - cand. sci. (techn.), senior researcher, Institute of Mechanics and Engineering, Kazan Science Center, Russian Academy of Sciences (IME KazSC RAS)

Yuriy A Volkov - cand. sci. (phys.-math.), general director, Oil field development methods improvement center ltd.

nocmynunu e pedu^uw 18 OKmnöpn 2017 z.