В работе представлены результаты исследований авторов по созданию пакета программ ТранSшельф, предназначенного для расчета транспортировки природного газа по морским трубопроводам с учетом специфики арктического шельфа и особенностей замерзания морской воды. Пакет позволяет рассчитать все характеристики потока и динамику оледенения (оттаивания) внешней поверхности трубопровода.
Ресурсы углеводородного сырья в Арктике, по оценке экспертов, составляют треть мировых запасов углеводородов [1]. Разработка газоконденсатных месторождений на арктическом шельфе России является важной стратегической задачей топливно-энергетического комплекса. Ее решение невозможно без развития отечественного научно-технического потенциала. Приведенные в настоящей работе математические модели и пакет программ ТранSшельф вносят вклад в разработку отечественных программных продуктов для решения задач освоения углеводородных ресурсов в Арктике.
Как отмечается в работе [2], при освоении арктических месторождений по ряду причин предпочтение часто отдается трубопроводной транспортировке добываемого сырья.
Для анализа проектных решений по созданию новых морских трубопроводов, для обеспечения безопасности и надежности их эксплуатации необходимо проведение комплекса научно-исследовательских работ, в который входит непростая задача создания адекватной математической модели подводного трубопровода. Эта модель (цифровой двойник) в режиме реального времени должна давать информацию о процессе транспортировки сырья и о работе оборудования. Знание зависимостей всех характеристик потока от изменяющихся во времени реальных условий позволяет обеспечить безопасность и надежность транспортировки добываемого сырья.
Широкое распространение зарубежных программных комплексов OLGA, PIPESIM, Fluent, Star-CD и других создает иллюзию того, что с их помощью можно решить любую практическую задачу. В действительности это не совсем так, поскольку математические модели и вычислительные алгоритмы, используемые в этих программах, не универсальны. Адекватность математической модели в реальных сложных многофакторных задачах является предметом самостоятельного исследования. Нетривиальным моментом является также выбор эффективного численного метода решения систем уравнений модели, обеспечивающего возможность расчета в режиме реального времени.
В нормативных документах Газпрома подчеркивается, что для каждой практической задачи проектирования и расчета газопровода требуется доработка и коммерческих пакетов, и лицензионных программ с учетом специфики задачи. В последние годы использование зарубежных программ осложнилось из-за введения санкций в отношении России.
В настоящей работе представлена часть результатов исследований авторов по созданию пакета программ ТранSшельф [3], позволяющего рассчитать как характеристики смеси газов, транспортируемой по морским трубопроводам, так и динамику оледенения (оттаивания) внешней поверхности трубопроводов в северных морях.
Математическая модель течения газа по трубопроводу
Наиболее содержательные одномерные математические модели нестационарного неизотермического течения природного газа по трубопроводу постоянного круглого сечения имеют общую структуру [3−7]. Они состоят из уравнения неразрывности, уравнения движения (баланса импульса), уравнения сохранения полной энергии, уравнения состояния газовой смеси и калорического уравнения. В уравнение движения входит коэффициент гидравлического сопротивления, являющийся функцией расхода газа, внутреннего радиуса трубопровода, коэффициента относительной шероховатости внутренней стенки трубопровода и коэффициента динамической вязкости газа. Уравнения движения и энергии содержат учет рельефа трассы прокладки трубопровода. В уравнение энергии входит осредненное по углу значение радиальной составляющей вектора потока тепла на внутренней поверхности трубопровода. Величина радиальной составляющей вектора потока тепла определяется из решения сопряженной задачи теплообмена между потоком газа и окружающей средой через многослойную боковую стенку трубопровода, на внешней поверхности которого возможно образование слоя льда. Таким образом, поток тепла на внутренней поверхности трубопровода зависит от плотности, толщины, коэффициента теплопроводности и удельной теплоемкости каждого слоя обшивки трубопровода и слоя образовавшегося льда, от температуры окружающей воды и условий обтекания трубопровода, определяющих величину коэффициента теплообмена наружной стенки трубопровода с окружающей водой.
Плотность, давление, температура, скорость, массовые плотности полной и внутренней энергий газовой смеси в одномерной модели являются функциями координаты z, направленной вдоль оси трубопровода, и времени t.
Допустимость одномерного описания процессов, при котором все характеристики являются осредненными по сечению трубопровода, исследовалась в нашей книге [8]. В ней приведены расчеты профиля скорости в турбулентном потоке сжимаемого газа для полуэмпирических моделей турбулентности Прандтля, Кармана и Новожилова–Павловского. Решение этой задачи позволило оценить погрешность замены двумерной модели течения газа одномерной моделью. Проведенное исследование привело к следующему выводу: при малых числах Маха, при малости относительных уклонов трассы, при больших числах Рейнольдса и при малых значениях коэффициента эквивалентной шероховатости допустим переход к одномерной модели течения. В задачах транспортировки газа по морским трубопроводам в северных морях обычно все перечисленные условия выполнены. Численные значения этих оценок представлены в книге [3].
Система уравнений математической модели течения газа по трубопроводу дополняется начальными и граничными условиями, соответствующими рассматриваемой задаче. Названная математическая модель и ее варианты исследованы в нашей книге [3], там же рассмотрены алгоритмы численного решения системы уравнений модели.
Адекватность математической модели в каждой практической задаче обеспечивается соответствующим выбором уравнения состояния транспортируемой смеси газов в заданном диапазоне изменений давления и температуры потока. Газовые компании корректируют коэффициенты уравнений состояния в соответствии с данными специально проведенных натурных экспериментов. Эти поправки нередко являются коммерческой тайной, охраняемой законодательством.
Адекватность математической модели зависит также от расчета коэффициента гидравлического сопротивления λ и суммарного коэффициента теплопередачи β. При наличии экспериментальных данных о давлении и температуре газа в ряде сечений трубопровода надежнее всего определять коэффициенты λ и β из решения задачи идентификации этих коэффициентов. Алгоритм решения задачи идентификации итерационным методом квазилинеаризации подробно описан в наших работах [9–11]. Он позволяет рассчитать как величины λ и β, так и оценить необходимую точность измерения температуры и давления газа на выходе из рассматриваемого участка. Методика процедуры идентификации обобщается на произвольное число участков, удовлетворяющих условию постоянства на них λ и β.
Алгоритм решения задачи идентификации коэффициентов λ и β оформлен в виде программного комплекса ТранSшельф-rev, входящего в пакет программ ТранSшельф [3].
Для замыкания математической модели течения газа по трубопроводу ее необходимо дополнить моделью теплообмена потока газа с окружающей средой.
Математическая модель теплообмена при отсутствии слоя льда
При сверхвысоких давлениях (более 20 МПа) для потока газа характерны большие числа Рейнольдса: Re ≈108. При этом интенсивность турбулентных пульсаций в потоке приводит к тому, что в радиальном направлении лимитирующей стадией теплообмена с внешней средой является теплопроводность через многослойную стенку трубопровода. В рамках одномерной модели теплообмен с окружающей средой учитывается интегрально.
В правую часть уравнения энергии входит слагаемое 
имеющее размерность мощности объемного источника (стока) внутренней энергии; это слагаемое выражается через R – внутренний радиус трубопровода и через q(z, t) – осредненную по углу радиальную составляющую вектора потока тепла на внутренней поверхности трубопровода в z-ом сечении в момент t. Тепловые условия на внешней поверхности трубопровода на расстояниях порядка R допустимо считать неизменными по z и t. Дополнительный пульсационный перенос внутренней энергии газа в направлении оси z пренебрежимо мал по сравнению с конвективным переносом внутренней энергии в этом направлении, поэтому допустимо считать, что q(z, t) зависит от z и t только параметрически через зависимости от z и t температуры газа и внешних условий. Для осесимметричных задач вплоть до момента возникновения льда величина q(z, t) определяется из решения системы нестационарных одномерных уравнений теплопроводности в областях многослойной боковой стенки трубопровода и в области эффективного теплового пограничного слоя воды при соответствующих граничных условиях сопряжения и начальных условиях. Оценка толщины δ* теплового пограничного слоя воды рассмотрена в книге [8]. В пределах теплового пограничного слоя передача тепла моделируется линейным уравнением теплопроводности. Величина δ* зависит от многих факторов, в частности от донных течений. Если трубопровод заглублен в грунт, для расчета среднего теплового потока от окружающей среды необходимо решать двухмерную задачу теплопроводности, в которой температура является функцией радиуса и угла [12, 13].
Отсутствие осевой симметрии приближенно можно учесть в рамках одномерной модели за счет выбора эффективного значения параметра δ*, при этом величину δ* надежнее всего определять по экспериментальным данным с помощью процедуры идентификации λ и β и последующего расчета δ* по найденной величине β [3].
Необходимым условием для начала оледенения внешней поверхности трубопровода является такое падение температуры газа в потоке, при котором температура на внешней поверхности трубопровода становится ниже температуры T⁕ фазового перехода «морская вода – лед». Однако выполнение только этого условия еще недостаточно для начала оледенения. Надо, чтобы поток тепла от окружающей воды был строго меньше, чем поток тепла от внешней поверхности трубопровода внутрь его. Это условие заведомо не выполняется в теплых морях, в которых температура воды значительно выше температуры фазового перехода. Иначе обстоит дело в северных морях, например, в некоторых районах Баренцева моря придонная температура воды близка к температуре фазового перехода, равной в соленой морской воде T⁕ = 271,236 К.
Математическая модель теплообмена при оледенении трубопровода
Уравнения теплопроводности в слоях обшивки трубопровода остаются теми же, что и в модели теплообмена при отсутствии льда. Начальным условием для температуры в каждом слое служит распределение температуры в этом слое в момент начала оледенения. Модель дополняется нестационарным одномерным уравнением теплопроводности в области слоя льда, имеющего изменяющуюся толщину y = y(t), условием неизменности температуры на границе фазового перехода и модифицированным уравнением Стефана. В отличие от пресной воды, в соленой морской воде при нарастании льда происходит отток солей из образующегося льда в слой воды, прилегающий к фронту замерзания. Кроме того, теплота плавления, плотность и теплоемкость морского льда, его коэффициент теплопроводности и температура T⁕ фазового перехода являются сложными функциями солености морского льда и его температуры [14, 15]. Все это приводит к искажению динамики оледенения, которая имела бы место в пресной воде. Суммарный эффект выражается в замедлении процесса оледенения. В рассматриваемой модели этот эффект учитывается введением в модифицированное уравнение Стефана дополнительного притока тепла к фронту оледенения, пропорционального скорости нарастания льда, коэффициент пропорциональности α, как и дополнительный приток тепла, в пресной воде равен нулю. Коэффициент α является эффективным параметром модели оледенения. Выбор α и выбор остальных теплофизических параметров нарастающего морского льда представляет собой самостоятельную непростую задачу, для решения которой кроме математической модели необходимо наличие экспериментальных данных. Методика выбора согласованного набора теплофизических параметров нарастающего морского льда и параметра α приведена в работе [16] и в книге [3]. Эффективность этой методики была подтверждена при расчетах скорости оледенения поверхностей в соленой воде по экспериментальным данным, предоставленным К.Е. Сазоновым, начальником лаборатории «Морской ледотехники» ФГУП «Крыловский государственный научный центр».
Модель оледенения и оттаивания внешней поверхности многослойного трубопровода в северных морях, алгоритм численного решения этой нестационарной задачи методом явного выделения фронта и результаты расчетов вариантов, представляющих практический интерес, приведены в работе [17]. Там же даны качественные и количественные оценки допустимости перехода к квазистационарному варианту модели оледенения (оттаивания) многослойных областей. Эти оценки имеют большое значение при создании эффективных вычислительных алгоритмов расчета неустановившихся режимов транспортировки газа по морским трубопроводам при их возможном оледенении.
Алгоритм численного решения задач оледенения поверхностей в соленой морской воде оформлен в виде программных комплексов ТранSшельф-Ste-1 и ТранSшельф-Ste-2, входящих в пакет программ ТранSшельф [3].
Общая модель неустановившегося неизотермического течения природного газа и оледенения морских трубопроводов в северных морях
Эта модель состоит из уравнения неразрывности, уравнения движения, уравнения энергии, уравнения состояния Редлиха–Квонга и следующего из него уравнения связи внутренней энергии с температурой и плотностью газа, полуэмпирического уравнения Коулбрука–Уайта для расчета коэффициента гидравлического сопротивления λ и одного из блоков расчета радиальной составляющей q(z, t) вектора потока тепла на внутренней поверхности трубопровода. Модель дополняется начальными и граничными условиями.
Проведенный анализ известных численных методов решения одномерных задач газовой динамики привел к выбору для численного решения системы уравнений модели модифицированной явной двухшаговой схемы Лакса–Вендроффа, которая в этих задачах оказалась предпочтительнее других численных схем, как по скорости счета, так и по простоте реализации. Этот алгоритм численного решения системы уравнений модели транспортировки природного газа и оледенения морских трубопроводов в соленой морской воде оформлен в виде программного комплекса ТранSшельф -t [3].
Для выбора параметров модели, теоретический расчет которых затруднителен, предлагается использовать программный комплекс ТранSшельф-rev [3].
В работе [18] и в книге [3] приведены результаты расчетов тестовых вариантов задачи о неустановившемся течении газовой смеси с преобладанием метана по морскому трубопроводу при температурах, характерных для Баренцева моря. В качестве начальных данных в этих расчетах использовались характеристики установившегося течения смеси газов, расчет которых, включая расчет слоя льда, приведен в работе [19]. Для одного из тестовых вариантов в работе [20] исследована возможность всплытия трубопровода в результате его частичного оледенения.
Выводы
Результатом проведенного исследования явилось создание вычислительного алгоритма общей модели транспортировки смеси газов по морским трубопроводам, который позволяет рассчитать для неустановившихся режимов все характеристики потока и динамику нарастания (таяния) морского льда на внешней поверхности трубопровода в северных морях.
Расчет включает учет рельефа трассы, условий вдоль трассы, состава газовой смеси, конструкции трубопровода, нестационарности процессов теплообмена.
На основе этого вычислительного алгоритма компанией ООО «Питерсофтваре» был создан пакет программ ТранSшельф для персональных компьютеров, состоящий из семи программ, каждая из которых реализует численное решение определенной задачи. Подробное описание пакета программ ТранSшельф приведено в книге [3].
Пакет ТранSшельф является отечественным программным продуктом, предназначенным для расчета транспортировки сырья по реальным трубопроводам с учетом специфики арктического шельфа и особенностей замерзания соленой морской воды. Он позволяет без обращения к коммерческим пакетам получать в режиме реального времени информацию обо всех характеристиках потока и о состоянии морского трубопровода.
Литература
1. Дмитриевский А.Н., Еремин Н.А., Шабалин Н.А., Кондратюк А.Т., Еремин А.Н. Состояние и перспективы традиционного и интеллектуального освоения углеводородных ресурсов арктического шельфа // Деловой журнал Neftegaz. RU. – 2017. − № 1 (61). – С. 32−41.
3. Г.И. Курбатова, Н.Н. Ермолаева, В.Б. Филиппов, К.Б. Филиппов. Проектирование газопроводов в северных морях. – Санкт-Петербург: Лань, – 2020. – 352 с.
4. Васильев О.Ф., Бондарев Э.А., Воеводин А.Ф., Каниболотский М.А. Неизотермическое течение газа в трубах. – Новосибирск: Наука СО, –1978. –128 с.
5. Селезнев В.Е., Клишин Г.С., Алешин В.В., Прялов С.Н., Киселев В.В., Бойченко А.Л., Мотлохов В.В. Численный анализ и оптимизация газодинамических режимов транспорта природного газа. – М.: УРСС, 2003. – 223 с.
6. Лурье М.В. Теоретические основы трубопроводного транспорта нефти, нефтепродуктов и газа. – М.: Недра, – 2017. – 478 с.
7. Курбатова Г.И., Филиппов Б.В., Филиппов В.Б. Неизотермическое турбулентное течение сжимаемого газа // Математическое моделирование. 2003. – Том 15. – № 3. – С. 92−108.
8. Г.И. Курбатова, Е.А. Попова, Б.В. Филиппов, В.Б. Филиппов, К.Б. Филиппов. Модели морских газопроводов. СПб.: С.-Петерб. гос. ун-т, – 2005. – 156 с.
9. Курбатова Г.И., Ермолаева Н.Н. Анализ чувствительности модели транспортировки газа по морским газопроводам к изменениям параметров модели // Вестник Санкт-Петербургского ун-та. – 2019. – Сер. 10. – Т. 14. – Вып. 1. – С. 47–61.
10. Ермолаева Н.Н., Курбатова Г.И. Параметрическая идентификация модели установившегося неизотермического течения газа по морскому газопроводу // Морские интеллектуальные технологии. 2017. – № 1 (35). – Том 1. – С. 8–13.
11. G.I. Kurbatova, N.N. Ermolaeva, V. Philipoff. Effective calculation methods of the gas flow characteristics and offshore gas pipeline glaciation // 2nd International Conference on Mathematical Methods & Computational Techniques in Science & Engineering. Murray Edwards College. – University of Cambridge: MMCTSE, 2018.
12. Galina I. Kurbatova, Vladimir A. Klemeshev. Effective numerical methods for calculating non-stationary heat and glaciation dynamic processes for offshore gas pipelines. // Energy, 2020, Vol. 205, 117995, https://doi.org/10.1016/j.energy.2020.117995.
13. Курбатова Г.И., Виноградова Е.М. Аналитические и численные решения двумерных задач теплопроводности и электронной оптики // Вестник Санкт-Петербургского ун-та. – 2021. – Сер. 10. – Вып. 4. – С. 345–352.
14. Доронин Ю.П., Хейсин Д.Е. Морской лед. – Л.: Гидрометеоиздат, – 1975. – 320 с.
15. Сазонов К.Е. Материаловедение. Свойства материалов. Методы испытаний.
Лед и снег. – СПб.: Изд-во РГГМУ, – 2007. – 195 с.
16. Ермолаева Н.Н., Курбатова Г.И. Нестационарная модель нарастания морского льда // Вестник Санкт-Петербургского университета технологии и дизайна. Сер. 1. Естественные и технические науки. – 2017. – № 1. – С. 3–8.
17. Курбатова Г.И., Ермолаева Н.Н., Никитчук Б.Я. Модели оледенения и оттаивания внешней поверхности морского газопровода в северных морях // Математическое моделирование. – 2019. – Том 31. – № 5. – С. 3–19.
18. Курбатова Г.И., Ермолаева Н.Н. Программные комплексы расчета параметров транспортировки газа и оледенения морских трубопроводов на арктическом шельфе // Научно-технический сборник ВЕСТИ ГАЗОВОЙ НАУКИ. 2019. – № 2 (39). – С. 80–90.
19. В.Б. Филиппов. Модель морского газопровода с учетом оледенения // Вестник Санкт-Петербургского ун-та. – 2004. – Сер. 1 – Вып. 3. – № 1. – С. 103–111.
20. Malkov, V., Kurbatova, G., Ermolaeva, N., Malkova, Y., Petrukhin, R. Analysis of the strength of sea gas pipelines of positive buoyancy conditioned by glaciation // AIP Conference Proceedings 1959(1):050019; https://doi.org/10.1063/1.5034647.
