Особенности математического моделирования систем добычи и транспорта природного газа в Арктической зоне России

ё Э.А.Бондарев, И.И.Рожин, К.К.Аргунова

Особенности математического моделирования систем добьни...

УДК 622.279.72:519.633:622.691.4

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

Э. А.БОНДАРЕВ, И.И.РОЖИН, К.К.АРГУНОВА

ФГБУН «Институт проблем нефти и газа Сибирского отделения РАН», Якутск, Россия

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

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

Как цитировать эту статью: Бондарев Э.А. Особенности математического моделирования систем добычи и транспорта природного газа в арктической зоне России / Э.А.Бондарев, И.И.Рожин, К.К.Аргунова // Записки Горного института. 2017. Т. 228. С. 705-716. DOI: 10.25515/РМ1.2017.6.705

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

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

Задача образования гидратов в газовой скважине и алгоритм ее численного решения.

Необходимый анализ выполнен в рамках математической модели образования гидратов при течении газа в трубах, предложенной в монографии [12] и модифицированной в последующих публикациях авторов [1, 3, 4, 8, 14]. В этой модели течение газа в трубе изменяющегося со временем проходного сечения описывается в квазистационарном приближении, так как скорость

ё Э.А.Бондарев, И.И.Рожин, К.К.Аргунова

Особенности математического моделирования систем добьни...

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

Для описания образования и отложения гидратов используется квазистационарная математическая модель [1, 3, 4, 8, 12, 14], в которой движение несовершенного газа в трубах описывается в рамках трубной гидравлики, а динамика образования гидрата - в рамках обобщенной задачи Стефана, в которой температура фазового перехода «газ - гидрат» существенно зависит от давления в потоке газа. В этой модели, основанной на законах сохранения массы и энергии для потока газа, уравнения неразрывности, движения и энергии газа сведены к системе двух обыкновенных нелинейных дифференциальных уравнений относительно давления и температуры:

ф . Л?М2 ,1Ч

Л=-р ^5,пф-;

Т -6 £ — т - т)--^

ах ах срМ ср

где р - давление; х - координата вдоль оси трубы; р^ - плотность газа; - - гравитационное ускорение; ф - угол наклона трубы, отсчитываемый от фиксированной горизонтальной плоскости; у - коэффициент гидравлического сопротивления; М - р ^ и 0 - постоянный массовый расход

газа; « - безразмерное поперечное сечение; «0 - размерное поперечное сечение до образования гидратов; Т - температура газа; Б - диаметр проходного сечения; а - суммарный коэффициент теплопередачи; ср - удельная теплоемкость газа при постоянном давлении; Те - температура окружающих горных пород; и - скорость течения газа.

Плотность газа и коэффициент дросселирования 8 связаны с давлением и температурой уравнениями

Р КТ2 (ъ

р- — —; 8 --1 — 1 ,

2ЯТ Срр \дТ )р

где 2 — 2(р, Т) - коэффициент несовершенства газа как эмпирическая функция, зависящая от отношения давления и температуры к их критическим значениям; Я — 8,314/ц - газовая постоян-

п

ная; ц — £ у/ ц / - молярная масса газовой смеси; у, и ^ - объемная доля и молекулярный вес /-го

/—1

компонента природного газа.

В работе [2] было показано, что при высоких значениях давления и температуры достаточно хорошо согласуется с экспериментальными данными уравнение Латонова - Гуревича [7]

р

Г г Л

2 —

0,173761п Т + 0,73 Т

Рс

+ 0,1-Р,

Рс

где критические давление и температура газовой смеси определяются по правилу Кэйа [15]:

п п

рс — £ у/рс/, Тс — £ У/ТС1, в которых рС1 и ТС1 - критические параметры /-го компонента природного газа.

/—1 /—1

Уравнение, описывающее изменение площади проходного сечения скважины « с течением времени, записывается в безразмерном виде:

Ж ЭЛБондарев, И.И.Рожин, К.К.Аргунова DOI: 10.25515/PMI.2017.6.705

v Особенности математического моделирования систем добьни...

d-=b2 ^Fttt -

от 1 - b2ln S

где т = 1 hTc/phqhDlt - безразмерное время; b2 = a2D0/41h;Th(p) = alnp + b - равновесная температура образования гидрата; b1 = а^0/41 h; t - время; A,h, ph - коэффициент теплопроводности и плотность гидрата; qh - удельная теплота гидратообразования; D0 - диаметр трубы до образования гидрата; а2 - коэффициент теплообмена между слоем гидрата и горной породой; a1 - коэффициент теплообмена между газом и слоем гидрата.

Эмпирические коэффициенты а и b находят путем аппроксимации кривой термодинамического равновесия гидратообразования, определяемой методом Слоана [16] по известному составу газа. В уравнении (4) все значения температуры отнесены к критической температуре газа Tc. Начальные условия для уравнений (1), (2) и (4) сформулируем в виде

p

В уравнении (4) коэффициент ai зависит от изменяющейся во времени площади проходного сечения трубы S. Для вывода соответствующей зависимости используется известная полуэмпирическая формула для коэффициента теплопередачи при турбулентном течении газа в трубах [11]

Nu = 0,023Pr°'43Re°'8,

где Nu = a1Dj 1 g - число Нуссельта; Pr = ngcp /1 g - число Прандтля; Re = uDpgjng - число Рей-

нольдса; % и Xg - динамическая вязкость и коэффициент теплопроводности газа.

Используя выражение для массового расхода газа и формулу (6), получим искомую зависимость в виде

^ = 0,023Pr043

1 g

f

M

,0,8

D 0П g

4 y,8 1

4 1 —.

S 0,9 v '

В тех сечениях скважины, где образуется гидратный слой, т.е. где безразмерная величина проходного сечения £ < 1, коэффициент теплообмена в уравнении (2) вычисляется по формуле (7). При этом температура горных пород Те заменяется на равновесную температуру гидратообразования Т%.

Уравнения (2) и (4) содержат температуру горных пород Те, определяемую из решения дифференциального уравнения теплопроводности, которое записывается в виде, удобном для численного решения методом сквозного счета, т.е. с коэффициентами теплопроводности и объемной теплоемкости грунта как кусочно-постоянными функциями температуры в окрестности фазового перехода «лед - вода»:

iST 1 d

Г ¿г \

C

dt г дг

V

дг

где С - объемная теплоемкость горных пород; г - радиальная координата; X - коэффициент теплопроводности; г1 - наружный радиус скважины; г2 - радиус теплового влияния.

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

дТ

X

дг

В области горных пород от забоя скважины до подошвы многолетней мерзлоты коэффициенты в уравнении (8) являются постоянными, и его решение может быть выполнено стандартными методами. В области многолетней мерзлоты эта задача осложняется, так как здесь необходимо учитывать фазовый переход «лед - вода». Для численного решения задачи типа Стефана используются методы, основанные на подходе, изложенном в монографии [13]. Для этого метода авторы работ [5, 10] разработали экономичную разностную схему сквозного счета со сглажива-

ё Э.А.Бондарев, И.И.Рожин, К.К.Аргунова

Особенности математического моделирования систем добьни...

нием разрывных коэффициентов в уравнении теплопроводности по температуре в окрестности фазового перехода. Эта схема характеризуется тем, что граница раздела фаз явно не выделяется, что позволяет использовать однородные разностные схемы. При этом скрытая теплота фазового перехода W = qphрг&г вводится как сосредоточенная теплоемкость в коэффициент C (Те). Здесь

^¡¡к - удельная теплота фазового перехода «лед - вода»; рг и аг - плотность и весовая влажность горных пород.

На условном радиусе теплового влияния принимаем условие тепловой изоляции

дТ

—^ — 0 при г — г2.

дг

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

Гг0 - Гх, 0 < х < L - Н;

Те — \

е [Тгг, L - Н < х < L, v 7

где Те0 - температура на забое скважины; Г - геотермический градиент; Ь - глубина скважины; Н - мощность многолетней мерзлоты; Т-- - температура мерзлых горных пород.

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

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

I. Задаются геометрические и физические параметры, а также начальные условия (5) и (11).

II. При фиксированной площади сечения по уравнениям (1)-(3) методом Рунге - Кутта 4-го порядка вычисляются давление р(х) и температура Т(х) газа в скважине.

III. Из уравнений (4) и (7), делая шаг по времени, находится новое значение площади проходного сечения. При этом координата х входит в уравнение (4) как параметр.

IV. Определяется температурное распределение в массиве горных пород, т.е. решается задача (8)-(11). Поскольку сглаженные коэффициенты в уравнении (8) зависят от температуры, получающаяся разностная задача будет нелинейной и ее решение находится методом простой итерации с использованием прогоночных алгоритмов.

На каждом временном шаге пункты П-ГУ повторяются.

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

Расчеты выполнялись при следующих значениях параметров, соответствующих двум месторождениям Республики Саха (Якутия):

1) Средне-Вилюйскому а — 5,82 Вт/(м2-К); Б0 — 0,1 м; ф — 90о; ? — 0,02; рк — 920 кг/м3; qh — 510000Дж/кг; Xк — 1,88 Вт/(м-К); X8 — 0,0307 Вт/(м-К); ср — 2300 Дж/(кг-К); Пе — 1,3 -10-5 Па-с;

Я — 449,4 Дж/(кг-К); р0 — 240 -105 Па; Т0 — 323 К; рс — 46,573 -105 Па; Тс — 205,239 К; а — 7,009 К; Ь — 178,28 К; £ — 2550 м; Н — 500 м; Те0 — 328 К; Г — 0,0277 К/м; Тг — 271,15 К; Трк — 273,15 К; qph — 334400 Дж/кг; состав газа (объемные доли), %: СН4 90,34, С2Н6 4,98, С3Н8 1,74, /С4Н10 0,22, пС4Н10 0,41, С5Н12 + 1,55, С02 0,28, ^ 0,48;

2) Отраднинскому Я — 438,3 Дж/(кг-К); Б0 — 0,146 м; р0 — 188,35 -105 Па; Т0 — 286,35 К; рс — 44,71-105 Па; Тс — 195,376 К; а — 6,635 К; Ь —182,951 К; £ — 2480 м; Н — 680 м; Те0 — 286,48 К; Г — 0,0085 К/м; состав газа (объемные доли), %: СН4 83,15, С2Н64,16, С3Н8 1,48, /С4Н10 0,17, пС4Н10 0,50, /С5Н12 0,12, пС5Н12 0,17, С6Н14 0,17, С7Н16+ 0,28, С02 0,07, N2 9,50, Н2 0,02, Не 0,21; остальные параметры имеют те же значения, что и в первом варианте.

ё Э.А.Бондарев, И.И.Рожин, К.К.Аргунова

Особенности математического моделирования систем добьни...

Видно, что при примерно равной глубине продуктивного горизонта состав природного газа, а также пластовые и геотермические условия этих месторождений существенно различны. Характеристики горных пород приняты одинаковыми, индекс (И соответствует талой зоне, индекс /- мерзлой зоне горных пород (см. таблицу).

На начальном этапе вычислялся оптимальный массовый расход газа, соответствующий минимуму тепловых потерь в отсутствие гидратного слоя. Для Средне-Вилюйского месторождения он оказался равным 9 кг/с, а для Отраднинского такой расход примерно соответствует предельно свободному дебиту скважины и поэтому расчеты проводились при массовом расходе 2,86 кг/с, что соответствует рабочему дебиту скважины 187000 м3/сут. В последующем вычислительном эксперименте варьировались начальные значения свободного сечения скважины и массовый расход.

Физические характеристики горных пород

Интервал глубин, м рг, кг/м3 Юг, дол. ед. Хи, Вт/

0-86 2000 0,120 1,69 1,93 2570 2310

86-Н 2000 0,120 1,62 1,86 2680 2420

Н-980 2300 0,060 2,00 - 2440 -

980-1831 2350 0,055 2,27 - 2420 -

1831-2561 2380 0,053 2,38 - 2420 -

2561-£ 2330 0,057 2,10 - 2440 -

Наиболее интересные результаты получены для Средне-Вилюйского месторождения. Они представлены на рис.1-5. Прежде всего отметим, что при оптимальном расходе гидратная пробка образуется вблизи устья скважины, а ее нижняя граница находится много выше подошвы многолетней мерзлоты, что хорошо видно на рис. 1, а, где точка пересечения температуры газа (кривые 2 и 3) и равновесной температуры гидратообразования (кривая 4) соответствует глубине 2550 - 2505 = 45 м, т.е. в этом интервале температура газа становится ниже температуры гидратобразования. Уменьшение проходного сечения сопровождается резким падением давления вблизи устья скважины (рис.1, б). При сопряженной постановке задачи значения температуры и давления газа в верхней части ствола скважины оказываются немного выше, а интервал образования гидратной пробки немного меньше (33 м), чем при постоянной во времени температуре горных пород (кривые 3 на рис.1, а и 1, б).

Рассмотрим изменение площади проходного сечения скважины для двух значений массового расхода, когда в начальный момент скважина свободна от гидратов, т.е. 5(0) = 1 (рис.2). Видно, что при сопряженной постановке время образования гидратных пробок существенно возрастает в отличие от случая (см. поверхности 2), когда температура окружающих горных пород считается неизменной (см. поверхности 1). С увеличением массового расхода возрастает длительность процесса полной закупорки скважины гидратами: для оптимального расхода она составляет примерно 423 ч при сопряженной и 251 ч при несопряженной постановке, а для меньшего расхода соответственно 13 и 9 ч. При этом наиболее интенсивное протаивание происходит вблизи подошвы мерзлоты, что объясняется сравнительно высокой температурой газа, а радиус протаи-вания примерно пропорционален массовому расходу: 1,3 м при 9 кг/с, 0,25 м при 2 кг/с (сравни рис.3, а и рис.3, б).

Ситуация становится менее предсказуемой, если в начальный момент скважина только наполовину свободна от гидратов, т.е. 5(0) = 0,5 (рис.4, 5). В этом случае оптимальный расход газа составляет 4,8 кг/с. При таком расходе гидратная пробка вблизи устья скважины образуется через 173 ч (рис.4, а, поверхность 2) для сопряженной задачи и через 109 ч (рис.4, а, поверхность 1) для несопряженной задачи. Нижняя граница пробки находится на глубинах 56 и 52 м соответственно для несопряженной и сопряженной задач. При этом в нижней части скважины от забоя до глубин 1415 м (для несопряженной постановки) и 1160 м (для сопряженной) площадь проходного сечения со временем возрастает, т.е. происходит разложение гидратов. Для несопряженной задачи в конце процесса отбора газа на уровне забоя свободным от гидратов остается сечение = 0,86 (см. поверхность 1 на рис.4, а). При меньшем расходе - от забоя до глубин 1428 м (для несопряженной задачи) и 1173 м (для сопряженной) площадь проходного сечения со временем также растет. Выше этих отметок толщина гидратного слоя со временем растет, образуя гидратную

ё Э.А.Бондарев, И.И.Рожин, К.К.Аргунова

Особенности математического моделирования систем добьни...

Т, К 320 300 280 260

2 1

4 4

\ 5 1

0

510

1020 1530 2040 х, м

б

р, МПа 20

15

10

0

2Ч\ 1

\

1

0

510

1020 1530 2040 х, м

Рис. 1. Изменение температуры (а) и давления (б) газа с глубиной для скважины Средне-Вилюйского

месторождения при М = 9 кг/с

1 - t = 0,34 ч; 2 - t = 251,3 ч (при несопряженной постановке); 3 - t = 422,9 ч (при сопряженной постановке); 4 - равновесная температура гидратообразования; 5 - начальная температура горных пород

а

5 1

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

0 4

0

б

1 .

0,8 .

0,6

0,4

0,2

0 0 -0

1275

2550

423

t, ч

2550

13,5

Рис.2. Изменение площади проходного сечения скважины Средне-Вилюйского месторождения по глубине и во времени при 5(0) = 1 и М = 9 кг/с (а); М = 2 кг/с (б) 1 - несопряженная постановка, 2 - сопряженная постановка

а

0

1785

б

м Грк м  1,3 0,28

0,97 0,21

0,65 0,14 0,32 0,07

0 0 0 1785

2295

2550 423

t, ч

2550

13,5

Рис.3. Динамика изменения радиуса протаивания вокруг скважины Средне-Вилюйского месторождения по глубине и во времени при 5

а

0

х, м

Рис.4. Изменение площади проходного сечения скважины Средне-Вилюйского месторождения по глубине и во времени при 5(0) = 0,5 и М = 4,8 кг/с (а); М = 2 кг/с (б) 1 - несопряженная постановка, 2 - сопряженная постановка

1785

2040

х, м 2295 2550 "173,2

б

м ГрИ, м 0,83 0,91

0,63 0,68 0,42 0,45 0,21 0,23

0 0 0 1785

Г, ч

2040

2295

X м 2550

69,7 139,4 , 209,1 Г, ч

Рис.5. Динамика изменения радиуса протаивания вокруг скважины Средне-Вилюйского месторождения по глубине и во времени при 5

а

0

Рис.6. Изменение площади проходного сечения (а) и радиуса протаивания (б) скважины Отраднинского месторождения по глубине и во времени при 5(0) = 1 и М = 2,86 кг/с

ё Э.А.Бондарев, И.И.Рожин, К.К.Аргунова

Особенности математического моделирования систем добьни...

пробку вблизи устья через 209 ч для сопряженной задачи и через 9 ч для несопряженной задачи (рис.4, б). Для несопряженной задачи в конце процесса отбора газа на уровне забоя свободным от гидратов остается сечение 5 = 0,53 (см. поверхность 1 на рис.4, б), а для сопряженной задачи -часть скважины от забоя до отметки 357 м полностью очищается от гидратов. Отметим, что эти отметки превышают глубины, на которых температура газа становится равной равновесной температуре гидратообразования. Глубина протаивания горных пород в этом случае меньше, чем для скважины со свободным от гидратов в начальный момент сечением, и для этих двух расходов не превышает 1 м (рис.5, а и 5, б). Это объясняется гораздо меньшим временем теплового воздействия газа на горные породы.

Теперь рассмотрим соответствующие процессы для Отраднинского месторождения. Это месторождение отличается от Средне-Вилюйского низкой пластовой температурой, которая близка к равновесной температуре гидратообразования. Из динамики изменения площади проходного сечения скважины (рис.6, а) следует, что если вначале проходное сечение скважины было свободно, то при расходе 2,86 кг/с полная закупорка гидратами происходит за 4,8 ч, при этом радиус протаивания горных пород составляет примерно 0,15 м (рис.6, б). Поверхности для несопряженной и сопряженной постановок задачи здесь практически сливаются, и поэтому эти обозначения не нанесены.

Моделирование работы магистрального газопровода. Рассмотрим этот процесс на примере магистрального газопровода «Сила Сибири», который предназначен для транспорта природного газа Чаяндинского и Ковыктинского месторождений, расположенных на территории Республики Саха (Якутия) и Иркутской области, в Дальневосточные районы России, а также -его экспорта в Китай и другие страны Азиатско-Тихоокеанского побережья. На первом этапе проектируется строительство участка «Чаянда - Ленск» протяженностью 207 км, диаметром труб 1,4 м с рабочим давлением 9,82 МПа. Этот участок характеризуется почти сплошным распространением многолетней мерзлоты. Для предотвращения нежелательных последствий ее протаи-вания в проекте предусмотрена теплоизоляция трубы из материала «Экстрол-45» (экструдиро-ванный пенополистирол) толщиной 0,2 м с коэффициентом теплопроводности 0,029 Вт/(м-К). Учитывая высокую стоимость материала (примерно 6900 руб./м3), целесообразно оценить необходимость такого проектного решения.

Дополнительно рассмотрим последствия возможного отклонения от проектного решения, которым предусматривается тщательная осушка газа перед подачей в магистраль. Дело в том, что пластовые условия Чаяндинского месторождения (температура 282,15 К, давление 13 МПа) [9] полностью соответствуют равновесным условиям образования гидратов природного газа данного состава, %: СН4 85,1366, С2Н6 4,5969, С3Н8 1,5641, /С4Н10 0,5886, /С5Н12 0,1734, С02 0,1441, К2 7,3031, Не 0,4034, Н2 0,0646, СН3ОН 0,0226. Его осушка, согласно проекту, должна довести мольную долю воды до 0,0026 %. Максимально добываемое количество газа на Чаяндинском месторождении должно составить 25 млрд м3/год, т.е. примерно 700 кг/с.

Для описания исследуемого процесса воспользуемся системой (1)-(2), в которой положим угол наклона трубы к горизонтальной плоскости равным нулю, а коэффициент несовершенства газа будем определять из уравнения Бертло в форме, предложенной в монографии [6]:

( т2

2 = 1 + 0,07^-

Рс Т

р Т. [. Т

1" 6-г 2

V Т J

Вычисления проводились при следующих исходных данных: 00 = 1,4 м, длина участка

L = 200000 м, рабочее давление р0 = 98105 Па, температура на входе Т0 = 282,15 К, коэффициент

теплообмена с учетом теплоизоляции а0 = 0,145 Вт/(м2-К), начальная температура грунтов Тр = 271,15 К; для слабольдистой супеси коэффициенты теплопроводности и объемные теплоемкости в талом и мерзлом состояниях соответственно X й = 1,6 Вт/(м-К), X р = 1,7 Вт/(м-К), Са = 2,8 -106 Дж/(м3-К), Сг = 2,1-106 Дж/(м3-К); плотность р = 1760 кг/м3, влажность ю = 0,233. Для природного газа Чаяндинского месторождения указанного состава найдены газовая постоянная R = 453,524 Дж/(кг-К), критические параметры рс = 45,01-105 Па и Тс = 195,075 К, коэффициенты

ё Э.А.Бондарев, И.И.Рожин, К.К.Аргунова

Особенности математического моделирования систем добьни...

а = 10,73 К и Ь = 117,979 К. Массовый расход газа М = 700 кг/с. Заглубление стальных труб марки К60 толщиной стенки 0,032 м и теплопроводностью 68 Вт/(м-К) предусматривается 1,5 м до верха трубы. Остальные параметры имеют те же значения, что в задаче образования гидратов в скважине.

а б

Рис.7. Распределение температуры (а) и давления (б) газа по длине полностью теплоизолированного

газопровода

1 - t = 4,2 мин; 2 - t = 70 ч; 3 - t = 120 ч; 4 - ТрН; 5 - Т„ при t = 120 ч а б

Рис.8. Распределение температуры (а) и давления (б) газа по длине частично теплоизолированного

газопровода

1 - t = 4,2 мин; 2 - t = 70 ч; 3 - t = 120 ч; 4 - ТрН; 5 - Т„ при t = 120 ч

Т, К 280 275 270 265

260

б

р, МПа 9 8 7

6

1

3 2 1

5104

Рис.9. Распределение температуры (а) и давления (б) газа по длине газопровода после 120 ч работы

для различных типов теплоизоляции 1 - частичная; 2 - по всей длине; 3 - без изоляции; 4 - Тр1,; 5 - Т„

а

1 105 1,5 105 х, м

5 104 1 105 1,5 105 х, м

0

0

ё Э.А.Бондарев, И.И.Рожин, К.К.Аргунова

Особенности математического моделирования систем добьни...

Реализовывались два сценария: 1) идеальное выполнение проекта, т.е. в газопровод подается осушенный газ, что полностью исключает образование гидратов; 2) в газопровод подается влажный газ, при этом наряду с другими параметрами (давление и температура газа, динамика оттаивания грунта) определяется динамика гидратообразования и расчет проводится до тех пор, пока давление на выходе не опустится ниже допустимого значения 4 МПа. В каждом из сценариев анализировались три варианта расчетов: 1) теплоизоляция трубы на всем участке; 2) теплоизоляция начального участка трубы, длина которого определялась по условию равенства температуры на контакте внешней поверхности трубы с грунтом температуре протаивания-промерзания грунтов Jph; 3) теплоизоляция отсутствует.

Анализ результатов вычислений начнем с первого сценария (рис.7-10). Основной вывод из представленных на этих рисунках кривых следующий: имеет место значительное охлаждение газа за счет дросселирования, что существенно снижает опасность протаивания мерзлого грунта. Более того, при полностью теплоизолированном трубопроводе газ охлаждается ниже температуры замерзания воды в грунте примерно через 100 км, но при этом температура на контакте газопровода с грунтом все время ниже этой температуры (рис.7, а), т.е. грунт не будет протаивать. Отметим также, что распределение температуры и давления в газопроводе очень быстро выходит на стационарный режим (кривые 1 на рис.7, а и 7, б) и в дальнейшем не изменяется, из-за чего кривые 1, 2, 3 на этих рисунках практически сливаются. Это объясняется резким различием длительности переходных процессов в газопроводе и в окружающем его грунте.

При частичной теплоизоляции (первые 90 км) контактная температура через 5 сут становится равной температуре фазового перехода Т^ только на головном участке (примерно 30 км), а затем существенно понижается (рис. 8, а). Отсюда следует, что затраты на теплоизоляцию можно сократить как минимум в 2 раза.

Как и следовало ожидать, без теплоизоляции температура газа на начальном участке немного ниже, чем при частичной теплоизоляции, а затем они сравниваются (см. кривые 1 и 3 на рис.9, а). В отсутствие теплоизоляции контактная температура превышает температуру фазового перехода Тр^ на начальном участке примерно в 25 км (кривая 5 на рис.9, а). В этом случае грунт оттаивает примерно на таком же участке, а радиус протаивания через 5 сут не превышает 0,8 м (рис.10). Как видно из рис.7, б - рис.9, б, температурный режим практически не влияет на распределение давления в газопроводе.

Перейдем к анализу результатов вычислений для второго сценария, в котором допускается, что по тем или иным причинам в трубопровод поступает влажный газ, так что при анализе параметров его транспорта следует учитывать возможность образования газовых гидратов и, следовательно, уменьшение проходного сечения трубы по мере нарастания гидратного слоя. Эти результаты представлены на рис.11-13. Вычисления продолжались до того времени, когда давление на выходе становилось равным 4 МПа.

Прежде всего отметим, что время вычислений при реализации второго сценария существенно сокращается, а это означает, что процесс образования гидратов проходит достаточно интенсивно, и за короткое время (около 7 ч) давление на выходе понижается до заданного значения. При этом температура газа всюду ниже равновесной температуры гидратообразования, и со временем это различие возрастает (сравни кривые 1 и 2 на рис.11, а и 12, а). Кроме того, в отличие от первого сценария давление со временем также понижается (сравни кривые 1 и 2 на рис.11, б и 12, б), что вызвано уменьшением проходного сечения трубы. Еще более важно, что в этом сценарии температура на контакте трубы с грунтом всегда ниже температуры протаивания грунта

х, м

Рис.10. Динамика протаивания грунта вокруг газопровода без теплоизоляции по его длине

ё Э.А.Бондарев, И.И.Рожин, К.К.Аргунова

Особенности математического моделирования систем добьни...

Рис. 11. Распределение температуры (а) и давления (б) газа по длине полностью теплоизолированного

газопровода

1 - t = 4,2 мин; 2 - t = 7 ч; 3 - Ти при t = 7 ч; 4 - ТрИ; 5 - Т„ при t = 7 ч

Т, К 290 280 270 260 250

4 5

2 1

1

0

5104

1 105 1,5 105 х, м

б

р, МПа 8 6 4 2

1

2

1

0

5104

1 105 1,5 105 х, м

Рис. 12. Распределение температуры (а) и давления (б) газа по длине газопровода без теплоизоляции 1 - t = 4,2 мин; 2 - t = 6,3 ч; 3 - Ти при t = 6,3 ч; 4 - ТрИ; 5 - Т„ при t = 6,3 ч

а

Тр, (сравни кривые 4 и 5 на рис.11, а и 12, а). Здесь проявляется эффект так называемой «тепловой завесы», роль которой выполняет образующий слой гидрата. Из-за этого эффекта распределение температуры по длине частично теплоизолированного газопровода и газопровода без изоляции оказывается практически одинаковым, и

поэтому здесь приводятся только результаты для последнего варианта расчетов (рис.12). Динамика распределения площади проходного сечения по длине трубы (рис.13) показывает, что влияние теплоизоляции на толщину гид-ратного слоя неоднозначно. При наличии теплоизоляции толщина этого слоя меньше, чем при ее отсутствии почти по всей длине (примерно 170 км) газопровода. Однако на конечном участке ситуация меняется на обратную (сравни поверхности 1 и 2 на рис.13).

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

5 1

0,97 -0,93 -0,9 -0,86 "

0,83 -0

шшш

110

х, м

1,510

2105

6,96

t, ч

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

0

ё Э.А.Бондарев, И.И.Рожин, К.К.Аргунова

Особенности математического моделирования систем добьни...

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

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

Результаты вычислительного эксперимента позволяют: 1) оценить риск аварийного снижения давления при недостаточной осушке газа; 2) снизить затраты на теплоизоляцию газопровода.

ЛИТЕРАТУРА

1. Аргунова К.К. Математические модели образования гидратов в газовых скважинах / К.К.Аргунова, Э.А.Бондарев, И.И.Рожин // Криосфера Земли. 2011. Т. 15. № 2. С. 65-69.

2. Аргунова К.К. Свойства реального газа и их аналитическое представление / К.К.Аргунова, Э.А.Бондарев, И.И.Рожин // Газохимия. 2010. № 6 (16). С. 52-54.

3.АргуноваК.К. Численное изучение нелинейных эффектов в моделях добычи природного газа: Автореф. ... канд. физ.-мат. наук / Якутский государственный университет им. М.К.Аммосова. Якутск, 2005. 18 с.

4. Бондарев Э.А. Математические модели образования гидратов в газовых скважинах / Э.А.Бондарев, К.К.Аргунова // Информационные и математические технологии в науке и управлении: Труды XIV Байкальской всероссийской конференции. Иркутск: ИСЭМ СО РАН, 2009. Ч. 3. С. 41-51.

5. Будак Б.М. Разностный метод со сглаживанием коэффициентов для решения задачи Стефана / Б.М.Будак, Е.Н.Соловьева, А.Б.Успенский // Журнал вычисл. математики и мат. физики. 1965. Т. 5. № 5. С. 828-840.

6. Вукалович М.П. Уравнение состояния реального газа / М.П.Вукалович, И.И.Новиков. М.-Л.: Госэнергоиздат, 1948. 340 с.

7. Латонов В.В. Расчет коэффициента сжимаемости природных газов / В.В.Латонов, Г.Р.Гуревич // Газовая промышленность. 1969. № 2. С. 7-9.

8. Определение интервала гидратообразования в скважинах, пробуренных в многолетнемерзлых породах / К.К.Аргунова, Э.А.Бондарев, В.Е.Николаев, И.И.Рожин // Электронный научный журнал «Нефтегазовое дело». 2008. http://www.ogbus.ru/authors/Argunova/Argunova_2.pdf. 11 с. (дата обращения 2.10.2017).

9. Перепеличенко В.Ф. Перспективы освоения уникального нефтегазоконденсатного месторождения Якутии // Электронный научный журнал «Георесурсы. Геоэнергетика. Геополитика». 2012. Вып. 1 (5). http://oilgasjournal.ru/vol_5/ perepelich.pdf. 8 с. (дата обращения 4.09.2017).

10. Самарский А.А. Экономичная схема сквозного счета для многомерных задач Стефана / А.А.Самарский, Б.Д.Моисеенко // Журнал вычисл. математики и мат. физики. 1965. Т. 5. № 5. С. 816-827.

11. Теория тепломассообмена / С.И.Исаев, И.А.Кожинов, В.И.Кофанов, А.И.Леонтьев, Б.М.Миронов, В.М.Никитин, Г.Б.Петражицкий, М.С.Самойлов, В.И.Хвостов, Е.В.Шишов. М.: Высшая школа, 1979. 495 с.

12. Термогидродинамика систем добычи и транспорта газа / Э.А.Бондарев, В.И.Васильев, А.Ф.Воеводин, Н.Н.Павлов, А.П.Шадрина. Новосибирск: Наука. Сибирское отделение, 1988. 272 с.

13. ТихоновА.Н. Уравнения математической физики / А.Н.Тихонов, А.А.Самарский. М.: Наука, 1977. 736 с.

14. Bondarev E.A. Modeling the formation of hydrates in gas wells in their thermal interaction with rocks / E.A.Bondarev, I.I.Rozhin, K.K.Argunova // Journal of Engineering Physics and Thermophysics. 2014. Vol. 87. N 4. Р. 900-907. DOI: 10.1007/s10891-014-1087-0.

15. Kay W.B. Density of hydrocarbon gases and vapors at high temperature and pressures // Industrial & Engineering Chemistry Research. 1936. Vol. 28. P. 1014-1019.

16. Sloan E.D. Clathrate hydrates of natural gases / E.D.Sloan, C.A.Koh. Boca Raton: Taylor & Francis Group/CRC Press, 2008. 720 p.

Авторы: Э.АБондарев, д-р техн. наук, главный научный сотрудник, bondarev@ipng.ysn.ru (ФГБУН «Институт проблем нефти и газа Сибирского отделения РАН», Якутск, Россия), И.И.Рожин, д-р техн. наук, ведущий научный сотрудник, i_rozhin@mail.ru (ФГБУН «Институт проблем нефти и газа Сибирского отделения РАН», Якутск, Россия), К.К.Аргунова, канд. физ.-мат. наук, старший научный сотрудник, akk@ipng.ysn.ru (ФГБУН««Институт проблем нефти и газа Сибирского отделения РАН», Якутск, Россия).

Статья принята к публикации 10.10.2017.