Проблемы учета профиля скорости в расчетах турбулентных течений в трубах Текст научной статьи по специальности «Математика»

УДК 532.517

Г. И. Курбатова, Е. А. Попова

Вестник СИбГУ. Сер. 10, 2007, вып. 2

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

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

Об одномерных и двумерных моделях транспортировки газа. К настоящему времени накоплен богатый опыт по расчетам транспортировки газа по трубам [6-9]. В начале этих исследований создавались модели, большинство которых было основано на раздельном описании гидравлических и тепловых процессов. Гидравлический расчет обычно проводился по одной из модификаций модели И. А. Чарного [7]. В упрощенном варианте одномерная нестационарная модель И. А. Чарного предполагает допустимость пренебрежения в уравнении движения всеми слагаемыми, кроме падения давления и трения. Однако столь упрощенное описание оказалось недостаточным для расчетов магистральных газопроводов высокого давления.

Важным этапом в моделировании транспортировки газа явилась одномерная нестационарная неизотермическая модель транспортировки неидеального сжимаемого газа по трубам, впервые предложенная и исследованная О. Ф. Васильевым, Э. А. Бондаревым, А. Ф. Воеводиным и М. А. Каниболотским [8] (далее модель ВБВК). В настоящее время она широко используется, например в работе В. И. Зубова, В. Н. Ко-терова, В. М. Кривцова и А. В. Шипилина [10], посвященной расчетам нестационарных газодинамических процессов в газопроводе на подводном переходе через Черное море. На этой модели основаны программно-математические комплексы «CorNet» и «AMADEUS» [11].

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

Наряду с одномерными рассматривались и двумерные модели неизотермических нестационарных и стационарных течений газа в трубах, учитывающие изменения в радиальном направлении, например [6]. Однако применительно к описанию транспортировки газа по морским трубопроводам при численном решении уравнений такого рода моделей возникает проблема, связанная с существенной разномасштабностыо процессов вдоль осей г и z.

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

Модель транспортировки газа при высоких давлениях. Приведем двумерную постановку задачи течения газа в трубах. Уравнения сохранения массы, баланса импульса и баланса внутренней энергии для установившегося осесимметричного неизотермического турбулентного течения по трубам многокомпонентной химически инертной смеси газов в терминах средних величин записываются в цилиндрической системе координат (r,0,z) следующим образом [1]:

© Г. И. Курбатова, Е. А. Попова, 2006

■уравнение неразрывности

^=0, (1)

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

ди dp 13,.

ри — = pgcosa(z) - — + - — (гг), (2)

¿/2 аг г or

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

де ди ди ,. _ ...

7Г = т "Т"" ~ Р1Г ~ dlv (3)

oz дг dz

Эта система дополняется приведенными ниже уравнением состояния р = р(р,Т) (см. (22)), калорическим уравнением е = е(р,Т) (см. (23)) и замыкающими уравнениями для касательной составляющей г тензора турбулентных напряжений (см. (4)) и для вектора плотности потока тепла q. В системе (1)-(3) плотность р, скорость и, внутренняя энергия е (соответственно температура Г), строго говоря, являются функциями двух переменных г и z. Только относительно давления р не вызывает сомнения возможность полагать р = p{z).

Рассмотрим возможность квазиодномерного подхода, при котором удается сохранить информацию о виде профиля скорости u(r,z) в каждом сечении, при том, что изменения всех термодинамических характеристик течения (р, Г, е, р) описываются на языке осредненных по сечению величин. Отметим, что такая возможность связана с физическими особенностями рассматриваемых задач. При сверхвысоких давлениях (около 23 МПа) средняя скорость транспортировки газа сравнительно невелика (примерно 3 м/с) и силы инерции малы, поэтому основной вклад в изменение плотности газа вносят термодинамические процессы, на которые большое влияние оказывают процессы теплообмена с окружаютцей средой. Заметные изменения температуры окружающей среды возможны только на масштабах порядка нескольких километров. Кроме того, заметное падение давления, вызванное трением, также наблюдается на этих масштабах, так как мала шероховатость внутренней поверхности используемых труб (для понижения шероховатости применяются специальные покрытия). Потому на масштабах трубы порядка сотен метров газ допустимо считать несжимаемым, что позволяет применять известные полуэмпирические модели турбулентности (Прандтля, Кармана, Новожилова-Павловского и др.) и, кроме того, пользоваться известным экспериментальным законом сопротивления (типа соотношения Коулбрука Уайта), который также получен для несжимаемых сред. На больших масштабах (порядка нескольких километров) величины р, Т, р претерпевают заметные изменения, которые представляют основной интерес для проектировщиков морских газопроводов.

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

Постановка задачи в квазиодномерном приближении. Один из подходов, казалось бы, мог состоять в пренебрежении зависимостью от г всех термодинамических величин (р, Т, е, р) и в сохранении зависимости от г только скорости и (и = u(r,z)). Однако такой подход применительно к системе уравнений (1)-(3) приводит к возможности лишь тривиального результата: Vz = о. На это обстоятельство впервые обратил внимание С. К. Матвеев [12]. Результат не столь очевидный, приведем его доказательство.

Известные полуэмпирические модели турбулентности Прандтля, Кармана, Новожилова-Павловского при условиях параметрической зависимости плотности от координаты 2 и при Vz = 0, как показано ниже (см. (18)), могут быть записаны в такой форме:

т = --Ф{(р(г)). (4)

Р

Здесь </?(r) = p(z)u(r,z) - интеграл уравнения неразрывности (1), вид функции ф зависит от вида полуэмпирической модели турбулентности. Существенно, что плотность в эти модели для г входит в виде множителя.

Утверждение. > Система уравнений (1), (2)

д(ри) = 0 dz

ди dp 13,.

ри~- = pg cos Ck(z) - — + гт),

oz az г or

т = НрЖФ)) (4')

при условии Mz = 0 приводит к равенству ^ = 0. <

Доказательство. Интеграл ри = </з(г) уравнения неразрывности при условии Mz ^ — 0 позволяет преобразовать конвективное слагаемое в уравнении движения (2) следующим образом:

ди 1 dp

Представление т (4') приводит к соотношению

-|-(гг)=/(р)~(н%>)). (6)

г or г аг

Подставим (5), (6) в уравнение движения (2) и разделим все члены на /(р), в результате получим

A(r)B(z) = D(z)+C(r), (7)

где A(r) = <p2(r), B(z) = -^d/z,

D(z) = - ^ +pgcosa(z)), C(r) = ~(гф(<р)).

Покажем, что в рассматриваемых задачах из уравнения (7) следует равенство А(г) = const. Действительно, существуют координаты Z\ ф z-2 такие, что B(z\) ф B(z-¿), так как в противном случае зависимость плотности от г носила бы фиксированный характер, что не имеет смысла в данных задачах. Из уравнения (7) следует система двух линейных уравнений относительно А{г) и С (г):

A(r)B(zl) = D(z1) + C(r), A(r)B(z2) = D{z2) + C{r).

Определитель системы отличен от нуля при B(z\) ф B(z-2). Решение существует и единственно, в частности для Л(г) находим

_ -D(zi) + D(z2) А{Г) - -B(Zl) + B(z2y (8)

Уравнение (8) с разделяющимися переменными, поэтому A(r) = const, т. е. ф2{г) = const. Отсюда следует, что Vz = 0 и u = u(z).

Вывод. Условие Vz = 0 совместно с системой (1), (2), (4) приводит к равенству = 0 Vz, которое позволяет записать осредненное по сечению одномерное уравнение движения в виде

dua , \ dp 2rw

pua-^r = pgcosa(z) - — + —, (9)

здесь tw(z) - значение касательного напряжения на внутренней стенке трубопровода, иа - средняя по сечению скорость потока, иа = и при ^ = 0 Уг.

Постановка задачи расчета функции iр(г). Упростим уравнение движения (2), воспользовавшись физическими особенностями задачи. Это позволит поставить задачу о приближенном расчете профиля ip(r) в сжимаемой среде. После того как она будет решена и функция ip(r) определена, перейдем к одномерному уравнению движения, осреднив уравнение (2) с учетом найденного профиля <р(г).

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

Характерная средняя скорость течения в исследуемых режимах, как отмечалось выше, мала и, следовательно, малы силы инерции. Уклоны трассы также незначительны. Эти обстоятельства позволяют заменить в уравнении движения (2) слагаемые /э«!^ и рд cos ct(z) на их средние по сечению значения, не внеся при этом большой погрешности. В результате получаем приближенное уравнение движения

9и ... dp Id,.

(ри— -рд cos a(z)) + — = - — (rr). (10)

R

Операция осреднения по сечению определена равенством (/) = f f(r)r dr.

о

Найдем левую часть в уравнении (10), осреднив для этого уравнение движения (2):

/ du ! \\ , dP z1 d f \\ 2Tw /1П'\

{ри _ _ рд COS a(z)) + - = (-— (rr)) = —. (10 )

С учетом (10') приходим к равенству

1(9 /А 2Г- ПП

которое совместно с замыкающим уравнением для г служит постановкой задачи об определении функции <р(г).

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

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

-Л

Р(и)2

(12)

зависимость ти, от г будет параметрической, обусловленной зависимостью от г величин ри {и). Коэффициент сопротивления А в данном подходе рассчитывается по одному из известных экспериментальных законов сопротивления, например по закону Коулбрука-Уайта

1

УА

= -216

к 2,51

774 шГТхУ'

(13)

В рассматриваемых задачах о течении сжимаемого газа интеграл уравнения неразрывности ц> = ри приводит к независимости от г не скорости, как в несжимаемых средах, а функции </9 = ри. Покажем, как в этом случае могут быть обобщены известные классические модели для т:

классическая модель Прандтля для несжимаемой среды

т = рГ-

в,и ¿г

йи , . .. г

/ = в1Я( 1-г), г=-;

обобщенная модель Прандтля

1

Р(*)

ФЛ<р(г))> ФЛ<р) = г

скр

¿г

скр ¿г '

модель Кармана для несжимаемой среды

т = р12

с1и

Шт

с1и дг'

I = ае2

йи / сРи

¿г / йг2

обобщенная модель Кармана 1

скр

скр ¿г

\с1г2 )

модель Новожилова-Павловского (модель Н~П) , ди

т = цгепТ£

дг

Т =

Р 1971

I | 2 >

I д г2 I

п = п(11е), эе„ = агп(11е);

обобщенная модель Н-П 1

Фз (<р(г)), Фз(<Р) =

ае„ (Ьр

уп-1 ¿г

й<р> з п , 2 я

¿г / &г2

Из равенств (14)—(16) следует, что обобщенные модели для т имеют вид

ф{<р(г)).

(14)

(15)

(16)

(17)

При выводе равенств (14)—(16) предполагается, что плотность является функцией только координаты г. Приближенное уравнение движения (11) совместно с равенствами (12), (17) приводит к следующей постановке задачи о расчете профиля функции

1-±{гфЫг))) = -^р2{и)\ (18)

Правая часть в (18) для установившихся течений - постоянная величина при постоянном Л. Действительно,

р{и) = (ри) =

здесь ф - массовый расход газа, Б = л Л2 - площадь поперечного сечения трубопровода.

Уравнение (18) - обыкновенное дифференциальное уравнение относительно ¡р(г), оно должно быть дополнено граничными условиями. В работах [2-4] найдены решения задачи о профиле функции <р(г) для различных обобщенных моделей т и при разных граничных условиях. В настоящей статье, пользуясь этими решениями, оценим влияние вида профиля скорости на характеристики течения, а также оценим погрешность замены профиля <р(г) на = сог^. Полученные ниже оценки справедливы в рамках сделанных предположений о малости сил инерции и уклонов трассы.

После того как задача расчета функции <р{г) решена, возможен переход от двумерной модели (1)-(3) к квазиодномерной постановке задачи расчета осредненных по сечению величин плотности, температуры, скорости, внутренней энергии. При этом удается сохранить, пусть и приближенно, информацию о виде профиля функции <р(г) и соответственно о виде профиля скорости. Осреднение по сечению уравнений движения (2) и баланса внутренней энергии (3) приводит к следующим уравнениям (их вывод приведен в [2]):

. с/ 1 Ар , . Лр(и)2 , ч

йТ _ Хр{и)2 р йр 2цш 5 3с йр

В (20) использованы уравнение состояния Редлиха-Квонга

_ НрТ___ср2

Р ~ 1 - 5р (1 + 5р)ТЧ* (21)

и согласованное с ним калорическое уравнение

йе _ (1Т Зс йр

2{1+5р)Т1/2 ~Иг' (22)

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

- значение радиальной составляющей вектора плотности потока тепла на внутренней поверхности газопровода при г = в г-ы сечении, расчет qw рассмотрен в работах [1, 2, 5].

Отличие квазиодномерного уравнения (19) от его одномерного варианта (9) заключается в записи инерционного слагаемого. В уравнении (9) инерционное слагаемое имеет вид

Аиа 2 2^1

^ар -г- = иар

аг аг р

и отличается от инерционного слагаемого в (19) тем, что (<р2) ф р2(и)2, по-

скольку (и2) Ф (и)2- Решение задачи о профиле функции ¡р(г) позволяет оценить относительную величину погрешности замены (</>2) на {ф)2 для различных режимов течения, удовлетворяющих условиям малости чисел Маха и относительных уклонов.

Приведем эти оценки для гидравлически гладких труб и режима развитой шероховатости.

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

Для режима гидравлически гладкого трубопровода в качестве модели т выберем обобщенную модель Н-П

1

, ае„

с/г

/>т'

Зга , ¿V •2 га

<1г / с/г2

с граничными условиями

I п ^

дт

— —оо,

(23)

(24)

г—И

п и ж - параметры модели Н-П. В режимах развитой шероховатости и переходном учет шероховатости в модели Н-П производится заменой граничного условия (24) на следующее условие [3]:

«и =

18

ди д?

= а» +К(кь, ¡у) 1/3,

(25)

г = 1

параметры; к - характерный линейный

ко-

где г = г/7? - безразмерный радиус; а*, Ь

размер шероховатости; и, - динамическая скорость = \/\тш\/р^; V = /х/р

эффициент кинематической вязкости.

Выбор параметров а», в граничном условии (25) для различных чисел Г1е достаточно громоздкий [3]. Как показали расчеты, в режимах развитой шероховатости и переходном целесообразнее использовать для г модификацию модели Кармана, основанную на изменении граничного условия в классической модели Кармана [2]. Обобщенная модифицированная модель Кармана

<1г

/¿V

^ (1т2

граничные условия

(26) (27)

здесь V, - динамическая скорость, аг-2 = 0,4 - постоянная Кармана, / - функция, содержащая учет шероховатости:

/ = - 23, 05

,11е /А

0,983

/0Л\ °>983

В =

7,75

17+ ЖеЛ/£

1.5 '

(28)

к - относительная шероховатость. 24

Для модели Н-П решение задачи (11), (23), (24) о профиле функции <р(г) следующее:

. / 2n — 1 \ (С\ [/ 2п-1\--гртг

4>(r,n) = [^j) [2) J dx> (29)

С =

J(n) = 1 j Beta(3z, у), у = (1 - п)/(1 + п), г = 2n/(2n - 1).

Здесь = "оРо; Ль «о - плотность и скорость газа при г = 0; n = n(Re) - параметр модели Н-П, его зависимость от числа Re найдена в работе [4] и приведена ниже:

Re 4 • 103 2,3 10" 1,1-105 1,1-10® 2 • 10® 3,2 • 106 Л 0,03950 0,02478 0,01747 0,01134 0,01027 0,00953 п 0,665 0,695 0,710 0,762 0,785 0,785

Значение коэффициента сопротивления Л определялось из закона сопротивления Коулбрука-Уайта (13).

Для обобщенной модифицированной модели Кармана решение задачи (11), (26)-(28) о профиле функции <р(г) имеет вид

АЛ, b-Vf „ /гЛ , , 1,25

ф) = 2,5<po\j\ (bln - (1 - >Гг)} , Ь = 1 + (30)

функция / определена равенством (28). Выражения для (ср2) и (<р)2 имеют вид

я к 2

(V2) = 4>2(г)г<1г, {<р)2 = ^у„(г)г*-) . (31)

о о

Для модели Н-П с учетом найденного решения (29) получаем

11 2 2 11 2 = л|г(У <*г) ¿г, (<р>2=2(1-х^)-1^ А^г) ,

(32)

О г

2 /2п- /СЛ^

Для характеристики отличия (<у32) от ((¿>)2 введем величину а - относительное отклонение (<р2) от (<р)2:

а = (<Р2) ~ Ы2

М2 '

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

0,01

0,7 0,75 0,8 п

Рис. 1. Зависимость относительного отклонения а от Не, рассчитанная по модели Н-П.

0,665 и достигается при Ие — 4 • 103. При этом числе Не величина а = 0, 04. Проведенные исследования вида профиля <р(г) при разных числах Рейнольдса и для различных коэффициентов шероховатости показывают, что величина а увеличивается при уменьшении 11е, так как профиль <р(г) при числах Ые < 105 существенно отличается от прямой.

Для модифицированной модели Кармана (</?2) и (ф)2 с учетом найденного решения (30) имеют вид

<у>2>=21 (33)

о

<<р)2 = 4 (I г ^2,1п ^^ - (1 - у^)) ) ¿гУ2. (34)

а

Рис. 2. Зависимость относительного отклонения а от Ие и к, рассчитанная по модифицированной модели Кармана.

Расчеты зависимости относительного отклонения а от числа Рейнольдса и относительной шероховатости для модифицированной модели Кармана представлены на рис. 2. Они показали, что в области развитой шероховатости максимальное значение а в исследуемом диапазоне изменений числа Рейнольдса и относительной шероховатости составляет 0,12.

При расчете интегралов, входящих в выражения (32), использовались специальные функции, а интегралы в (33), (34) берутся аналитически.

Заключение. Профиль функции ip (соответственно профиль скорости) входит в общую модель процессов (19)—(22) через величину {р2), являющуюся множителем в инерционном слагаемом уравнения движения. Рассчитанные значения отклонения (ip2) от (ip)2 в широком диапазоне изменений чисел Re и относительной шероховатости позволили оценить погрешность перехода от квазиодномерного к одномерному описанию рассматриваемых процессов.

Как показано на рис. 1 и 2, относительное отклонение п является малым в исследуемых режимах течения, поэтому следует ожидать, что с достаточной точностью можно пользоваться более простой одномерной моделью. Проведенные расчеты по квазиодномерной модели (11), (19), (20) и по одномерной модели (9), (20) подтвердили данное положение и позволили сделать

вывод о допустимости перехода к одномерному описанию исследуемых процессов транспортировки газа по морским газопроводам.

Summary

Kurbatova G. I., Popova E. A. Problems of velocity profile account for gas pipeline flow modeling.

Quasi-1D model approximately taking into account the velocity profile influence in turbulent gas flows in pipelines is presented. The error due to the transition to ID description in a wide range of Re and the relative roughness is calculated.

Литература

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

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

3. Филиппов Б. В., Филиппов В. Б. Профиль скорости турбулентного течения сжимаемого газа в шероховатых трубах// Нелинейные проблемы механики и физики деформируемого твердого тела/ Под ред. К. Ф. Черныха. СПб.: Изд-во С.-Петерб. ун-та, 2001. С. 73-84.

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

5. Филиппов В. Б. Расчет транспортировки газа по морскому газопроводу с учетом оледенения// Физическая механика/ Под ред. Б. В. Филиппова. СПб.: Изд-во С.-Петерб. ун-та, 2004. Вып. 8. С. 45-62.

6. Грачев В. В., Щербаков С. Г., Яковлев Е. И. Динамика трубопроводных систем. М.: Наука, 1987. 439 с.

7. Чарный И. А. Неустановившееся движение реальной жидкости в трубах. М.: Недра, 1975. 296 с.

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

9. Справочник по проектированию магистральных трубопроводов / Под ред. А. К. Дер-цакяна. Л.: Наука, 1977. 519 с.

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

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

12. Кочерыженков Г. В., Матвеев С. К. Приближенные методы расчета течений в трубах, каналах и щелях// Аэродинамика/ Под ред. Р. Н. Мирошина. СПб.: ВВМ, 2005. С. 134-151.

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