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

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

Сер. 10. 2010. Вып. 1

УДК 532.517

Е. В. Груничева, Г. И. Курбатова, Е. А. Попова

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ НЕСТАЦИОНАРНОГО НЕИЗОТЕРМИЧЕСКОГО ТЕЧЕНИЯ СМЕСИ ГАЗОВ ПО МОРСКИМ ГАЗОПРОВОДАМ

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

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

1) колебания газопотребления (суточные, сезонные и др.);

2) перевод системы на новый установившийся режим работы;

3) процессы пуска и перекрытия газопроводов;

4) обнаружение мест несанкционированного отбора и утечки газа;

5) оледенение подводной части газопровода в северных морях;

6) аварийные ситуации, сопровождающиеся разрывом линейной части газопровода.

Оперативное управление работой газопровода сводится к управлению неустано-

вившимися процессами. Ограничим исследование медленно меняющимися процессами и рассмотрим безударное неустановившееся неизотермическое турбулентное течение

Груничева Екатерина Викторовна — аспирант кафедры вычислительных методов механики деформируемого тела факультета прикладной математики-процессов управления Санкт-Петербургского государственного университета. Научный руководитель: проф. Г. И. Курбатова. Количество опубликованных работ: 4. Научное направление: математическое моделирование гидродинамических процессов. E-mail: mag_katris@rambler.ru.

Курбатова Галина Ибрагимовна — доктор физико-математических наук, профессор кафедры вычислительных методов механики деформируемого тела факультета прикладной математики-процессов управления Санкт-Петербургского государственного университета. Количество опубликованных работ: порядка 70. Научные направления: математическое моделирование, приложения тензорной алгебры, турбулентные течения. E-mail: gi_kurb@mail.ru.

Попова Елена Анатольевна — кандидат физико-математических наук, ассистент кафедры высшей математики факультета фундаментальных и гуманитарных дисциплин Санкт-Петербургского государственного горного института. Количество опубликованных работ: 13. Научное направление: математическое моделирование гидродинамических процессов. E-mail: e-lena-05@mail.ru.

© Е. В. Груничева, Г. И. Курбатова, Е. А. Попова, 2010

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

В работах [1, 3] построена квазиодномерная модель течения газовой смеси по морским газопроводам в установившемся режиме, учитывающая все перечисленные выше особенности задачи. Модель использовалась при проектировании морских газопроводов в Балтийском и Баренцевом морях (в частности, Северо-Европейского газопровода в Балтийском море и газопровода от Штокмановского газоконденсатного месторождения до Териберки в Баренцевом море). Предложенный в этих работах квазиодномерный подход позволил учесть влияние профиля скорости на основные характеристики течения. В [5] получено решение задачи расчета профиля скорости течения для полуэмпири-ческой модели турбулентности Новожилова-Павловского и модифицированной полуэм-пирической модели турбулентности Кармана, а также проведено сравнение результатов расчетов по одномерной и квазиодномерной моделям течения для установившихся режимов. В широком диапазоне изменений значений числа Рейнольдса и относительной шероховатости рассчитана величина относительной погрешности, возникающей при переходе от квазиодномерной к одномерной модели, и доказано, что при уменьшении числа Рейнольдса величина этой погрешности растет. Обоснован вывод о допустимости одномерного описания процессов транспортировки газа по морским газопроводам при условии, что характерное число Рейнольдса превышает величину 107.

К настоящему времени накоплен богатый опыт по расчетам транспортировки газа по магистральным газопроводам [6, 7]. Существенный вклад в решение этих задач внесли И. А. Чарный, Н. Н. Яненко, О. Ф. Васильев, С. К. Годунов, Э. А. Бондарев, А. Ф. Воеводин, М. А. Каниболотский и многие другие исследователи. Одномерная нестационарная модель транспортировки неидеального сжимаемого газа по трубам, впервые предложенная в работе [6], широко используется и в настоящее время, например, на ней основаны модель работы [7], посвященной изучению нестационарных газодинамических процессов в газопроводе на подводном переходе через Черное море, и программно-математические комплексы <Со^е1> и <AMADEUS> [8]. Проведенные нами исследования, как отмечалось выше, позволили оценить границы допустимости одномерного подхода. Предложенная ниже нестационарная математическая модель отличается от модели Васильева, Бондарева, Воеводина и Каниболотского [6] уравнением состояния и калорическим уравнением. Это связано с необходимостью возможно более точного описания термодинамических процессов при сверхвысоких давлениях.

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

Уравнение неразрывности

др дри

уравнение движения

дри д , 2 ч Хри2

~яГ + + р) =—-гб~ + Р9 сов а(г); (2)

дЬдг ' 4Е

уравнение Коулбрука-Уайта

1 „( к 2,51 \

7а~“ Чм +( }

43

Ие

2дс

пЕрс 1

уравнение баланса полной энергии

д(ре) д

дЬ дг

ри

Х° Г)

д, д0

пЕ рсис;

Е2Ш

(Т*(г, £) — Т)+ рид сов а(г),

(4)

3=1 3

+Г1п 1 + -

Гз-\) Х* \ гт

уравнение состояния Редлиха-Квонга

3—1

г 3 — 1 Е + ^ ^ ;

9=1

(5)

р =

ср

НрТ

1 -5р ~ (1 + 5р)Т1/2'

(6)

н = Ед/М, М = 2_^ пг т 1

к

Е

1

5 = Ь/М, с = а/М2, Ь = ОъЕд Тс/рс,

Пг = 1,

Па(Ед )2Т2’5/рс.

(7)

В системе (1)-(7) приняты следующие обозначения: и, р, р, Т - скорость, плотность, давление и температура газа соответственно, являющиеся функциями времени £ и координаты г, направленной вдоль оси трубопровода; е = е(г,Ь), е = е(г,Ь) - массовые плотности внутренней и полной энергии соответственно; Е - внутренний радиус трубопровода; Хз, - коэффициент теплопроводности и толщина j-го слоя обшивки

(они могут быть функциями г, т. е. на разных участках конструкция обшивки может быть различной); т - количество слоев обшивки трубы; Х*, й* - коэффициент теплопроводности и эффективная толщина слоя льда на внешней поверхности газопровода; д - ускорение силы тяжести; а(г) - угол между направлением силы тяжести и осью трубопровода в г-м сечении; Х - коэффициент гидравлического сопротивления (в общем случае он может параметрически зависеть от г и £ через зависимости ^с = ^с(г, £) и рс = рс(г,1)); Ие - характерное число Рейнольдса; Qс, рс(г,1) - характерное значение массового расхода и динамического коэффициента вязкости газа соответственно; ка - коэффициент эквивалентной равномерно-зернистой шероховатости внутренней поверхности газопровода; Т* (г,£) - температура окружающей воды; Ед - газовая постоянная; тг, Пг - молекулярный вес и доля г-й составляющей газовой смеси соответственно; к - количество компонент газовой смеси; Оа, Оъ - числа, определяемые по значениям критических температуры Тс и давления рс согласно таблицам, приведенным в [9], для заданного химического состава газовой смеси.

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

Система уравнений (1)-(7) должна быть дополнена калорическим уравнением, задающим связь внутренней энергии газа с его температурой и плотностью. В данной работе внутренняя энергия е принята локально функцией независимых термодинамических переменных (Т, У), где V = ^ - удельный объем газа. Это позволяет записать равенство

де де де

Ле = + ду^ = С^Т + ду^ (8)

к

2

к

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

Су Су (Т).

В работе [3] получено выражение для производной

д^ = _ ,Т&Р дУ Р дТ’

которое для уравнения состояния Редлиха-Квонга приводит к равенству

де 3 ср2

дУ 2(1+6р)л/Т'

(9)

В движущейся среде равенство (8) задает выражение для материальной производной от внутренней энергии

д!е _ де д!Т де сГУ

скЬ дТ скЬ дУ скЬ ’

где ^ = -щ + и-т^. В силу уравнения неразрывности (1), материальная производная равна

^У = ^_=1ди (10)

<М скЬ р дг

Совместно с равенством (9) выражение (10) определяет выражение для материальной производной от внутренней энергии

(1'е д!Т 3 ср ди

Л Су <М +2(1 + 6р)у/т дг' 11

Равенство (8) позволяет выразить внутреннюю энергию е как функцию плотности и температуры газа с точностью до аддитивной постоянной ео:

т V

де 3 С

£ = J + ! о^^У = суТ — - 1п(1 + 5р) + ео- (12)

То Vo

Уравнение энергии (4) инвариантно относительно замены е на е + ео, т. е. не зависит от выбора начала отсчета внутренней энергии.

В системе уравнений (1)-(7) независимыми искомыми функциями являются р, и, е, р. Правая часть уравнения энергии (4) записана в терминах температуры, это же относится и к выражению (6) для давления. Поэтому необходимо найти зависимость температуры Т от плотности р и внутренней энергии е. Для этого воспользуемся уравнением (12), приведя его предварительно к безразмерному виду и сохранив для безразмерных величин прежние обозначения. Приведем искомую зависимость, опустив промежуточные выкладки:

Т(р,е) = ^со82 ^агссов (^зУ)) • (13)

Выражения для безразмерных комплексов и т* приведены ниже (см.(23)).

В ряде задач целесообразно вместо модели I использовать эквивалентную ей модель II, отличающуюся заменой уравнения для полной энергии (4) на эквивалентное ему уравнение баланса внутренней энергии. Получим его, вычитая из уравнения для полной энергии (4) уравнение движения (2), домноженное на и. В результате ряда преобразований приходим к уравнению

d>£ , ди 2 №*/ ^ ал , А/ж3

^+р^ = я^(Г (M)-T) + 7ir

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

d'f _ dpf_ dpuf

^ dt дt dz ’

в результате получим тепловое уравнение, которое, в силу остальных уравнений системы (1)—(7), эквивалентно уравнению полной энергии (4):

(дрТ дриТ\ ди 3 ср2 ди 2 Хри3

с'{-Ж + ^г) = -ра; - 2(1 + ад^ + mv(T (zJ) - Г) + Т7Г (14)

Система уравнений (1)—(7) дополняется начальными и граничными условиями, соответствующими исследуемой нестационарной задаче.

В качестве примера приведем условие для задачи о выходе на новый установившийся режим при заданном изменении расхода на правом конце газопровода. Начальными данными в ней служат параметры установившегося режима po(z), To(z),

t = 0: pu = const, p(z, 0) = p0(z), T(z, 0) = T0(z). (15)

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

В данной задаче на левом конце газ втекает со скоростью, меньшей скорости звука, при этом, как известно, число условий должно быть равно m +1, где m - размерность задачи. На правом конце газ вытекает со скоростью, меньшей скорости звука, при этом ставится одно граничное условие. В рассматриваемой задаче граничные условия имеют вид

t> 0,z = 0: p(t, 0) = const = p0, T (t, 0) = const = T0, ( )

t> 0,z = L: pu = y = y* (t), ( )

здесь p0 и T0 - заданные неизменные плотность и температура поступающей на вход смеси газов соответственно; y*(t) - известный закон изменения удельного расхода газа на правом конце газопровода.

Приведем уравнения моделей I и II, а также начальные данные (15) и граничные условия (16) к безразмерному виду. Введем характерные масштабы длины zx, времени tx, скорости ux, плотности px, давления px, температуры Tx и внутренней энергии ех. Перейдем к безразмерным переменным и функциям, сохранив для соответствующих безразмерных величин прежние обозначения. После преобразований уравнения моделей I и II, дополненные (15) и (16) в безразмерной форме, запишутся следующим образом:

Безразмерная форма модели I

dp dpu

i + ~к = °- (17)

dpu д , , \ 2

-7^-+ (yOM + rriip) = -т2ри +тзрсова{г), (18)

д f pu2\ д ( pu3 ^

— l,oe + т10— I + — I рие + т10— + т5ир ]

= m6 (T*(z, t) — T(p,e)) + mnpu cos a(z), (19)

PT(P’£) P2

p = ”'FM'"'*(i + f.(,jr'%£)’ (20)

t = 0: pu =1, p = p0(z), T = T0(z),

t> 0,z = 0 : p = 1, T =1, (21)

t > 0,z = L : pu = y*(t).

Зависимость T(p,e) определена равенством (13).

Безразмерная форма модели II

Соотношения (17), (18), (20), (21) остаются без изменения, а уравнение (19) заменяется безразмерной формой уравнения баланса тепла (14)

д(pT) d(puT) f m7p2 \ du

» + ~аГ~ + \’щр+ (1 + адг-/з) Э1 = - Т) + ,щ (22)

Безразмерные комплексы т\ - тц, т*, 6* в уравнениях (18)—(20), (22) выражаются через характерные значения и параметры задачи по формулам

рх Хгх гх д Хи2, гх

'-'x'J ''^x^x

'т-\ = --------^ , то = —— , точ = —, тл = ------------------------------------------------—— ,

рхи2х 4 R м2 4 cvRTx

px 2zx 3cpx hpxTx

m5 = -----------7fT , m6 = —^-------------------------------------------------------------— , m7 = --r-- , m8 —

pxCvTx’ cvR2pxuxW ’ 2cvT3/2 ’ Px

22 cp2x ™ _ux _ _zx9 x _ г _ Px

гх х хс> г* г*

Ш9 = -----> ШЮ = -------- ? ™11 =----- , Ь* = ЬрХ, — , •

РхТх' ех ех рхех6*

Здесь Ь х = хх1их, ех = суТх, безразмерный комплекс т* входит в уравнение (13).

Вычислительные алгоритмы. Системы уравнений моделей I и II, записанные в эйлеровых координатах, решались численно методом конечных разностей. Вводилась равномерная сетка ^+1 = Х1 + Дх и £п+1 = Ьп + Д(, где Дх и ^ - величины шагов по координате г и по времени соответственно. Рассматривались различные варианты аппроксимации дифференциальных уравнений конечно-разностными соотношениями. В основном применялись двухслойные неявные разностные схемы, имеющие второй порядок аппроксимации по пространству и первый - по времени. Система разностных уравнений - замкнутая нелинейная алгебраическая система уравнений относительно искомых сеточных функций на новом временном слое. Решение строилось с помощью

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

1

где норма \\1 я\\ определялась как максимум модуля по всем значениям /*, г = 0, 1,..., N — 1; N - число точек сетки по координате г; в - номер итерации. Начальными данными служили параметры установившегося режима, рассчитываемые по соответствующей системе дифференциальных уравнений методом Рунге-Кутта.

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

Для всех вычислительных алгоритмов оптимальная величина шага по времени А* определялась в результате численного эксперимента. Так, в первом алгоритме при оптимальном значении шага А* требовалось не более трех итераций во всех итерационных процессах при условии ерв = 10~7.

Проведенные расчеты многочисленных вариантов задачи о выходе на новый установившийся режим подтвердили известное положение о преимуществе дивергентной формы записи уравнений системы при численном решении. Так, при использовании модели II, независимо от вычислительного алгоритма, для некоторых вариантов задания граничной функции у*(Ь) возникали искажения решения нефизического характера, которые не устранялись измельчением сетки. Но при переходе к модели I эти проблемы снимались.

Проведенные исследования позволяют сделать следующий вывод: параметры установившегося режима ро(г), То(г), ро(г) можно рассчитывать численно по стационарному варианту модели II; численное решение нестационарных задач в общем случае должно базироваться на модели I.

Литература

1. Курбатова Г. И., Филиппов Б. В., Филиппов В. Б. Неизотермическое турбулентное течение сжимаемого газа // Математическое моделирование. 2003. Т. 15, № 3. С. 92—108.

2. Дерцакян А. К., Курбатова Г. И., Неизвестнов Я. В., Филиппов Б. В. Некоторые научнотехнические проблемы освоения шельфа арктических морей России // Труды XIII сессии Междунар. школы по моделям механики сплошных сред. СПб., 1996. С. 99—109.

3. Курбатова Г. И., Попова Е. А., Филиппов Б. В. и др. Модели морских газопроводов. СПб.: С.-Петерб. гос. ун-т, 2005. 156 с.

4. Курбатова Г. И., Павловский В. А., Попова Е. А., Филиппов В. Б. О замыкающих уравнениях в моделях установившихся турбулентных течений в трубах // Вестн. С.-Петерб. ун-та. Сер. 1: Математика, механика, астрономия. 2003. Вып. 4 (№ 19). С. 76—88.

5. Курбатова Г. И., Попова Е. А. Проблемы учета профиля скорости в расчетах турбулентных течений в трубах // Вестн. С.-Петерб. ун-та. Сер. 10: Прикладная математика, информатика, процессы управления. 2007. Вып. 2. С. 18-28.

6. Васильев О. Ф., Бондарев Э. А., Воеводин А. Ф., Каниболотский М. А. Неизотермическое течение газа в трубах. Новосибирск: Наука, Сиб. отд., 1978. 128 с.

а

7. Зубов В. И., Котеров В. Н., Кривцов В. М., Шипилин А. В. Нестационарные газодинамические процессы в газопроводе на подводном переходе через Черное море // Математическое моделирование. 2001. Т. 13, № 4. С. 58-70.

8. Селезнев В. Е., Клишин Г. С., Алешин В. В. и др. Численный анализ и оптимизация газодинамических режимов транспорта природного газа. М.: УРСС, 2003. 223 с.

9. Рид Р., Праусниц Дж., Шервуд Т. Свойства газов и жидкостей / пер. с англ.; под ред. Б. И. Соколова. Л.: Химия, 1982. 532 с.

Статья рекомендована к печати проф. Ю. М. Далем.

Статья принята к печати 24 сентября 2009 г.