CC BY
32
3. Bondarchuk S.S., Vorozhtsov A.B. Mathematical Simulation of a Large Size Rocket Motors // Proceedings of the 3-rd International High Energy Materials Conference and Exhibit HEMCE-2000 (6-8 December 2000). Thiruvananthapuram, India.
4. Bondarchuk S.S. et al. Ignition Pressure Transient in Solid Rockets Initially Filled with Water // Journal of Propulsion and Power. V. 15. 1999. №6.
5. Bondarchuk S.S. et al. Solid Propellant Underwater Ignition Modelling // AIM Paper 99-0864 (A99-16700), AIAA, Aerospace Sciences Meeting and Exhibit, 37m, Reno, NV (11-14 Jan. 1999).
6. Bondarchuk S.S. et al. Methods for Recovery of Ingredients from Solid-Propellant Rocket Motors // The First International Symposium on Special Topics in Chemical Propulsion and Combustion of Energetic Materials, 18-22 June, Stresa (Lake Maggiore), Italy.
7. Godunov S.K. (ed) Numerical Solution of Multidimensional Problems of Gas Dynamics. M., 1982.
8. Vilyunov V.N., Zarko V.E. Ignition of Solids. Elsevier, Amsterdam, 1989.
9. Bondarchuk S.S., Vorozhtsov A.B. Gasdynamics of filling of elastic stagnation zones of large-sized solid rocket motors // 30*1 International Annual Conference of ICT (Energetic Materials) (29 June - 2 July 1999). Fraunhafer - Institute Chemische Technologie. Karlsruhe, Fed. Rep. of Germany.
10. Bondarchuk S., Vorozhtsov A., Salko A., Kondratova O. Mathematical Simulation of Airbag Inflation by Low Temperature Gas Generator Products // An International Journal Propellants. Explosives, Pyrotechnics. V. 25, Issue 5, 2000. (Published Online: 3 Nov. 2000 WILEY-VCH Verlag GmbH, Weinheim, Fed. Rep. of Germany).
11. Vorozhtsov A.B., Bondarchuk S.S., Salko A.E., Zima V.P. Inflation of Elastic Airbag under Wather - Experiment and Calculations // Energetic Materials: Ignition, Combustion and Detonation. 32nd International Annual Conference of ICT. (3 July - 6 July 2001). Karlsruhe, Fed. Rep. of Germany.
УДК (519.9+518.5):532.54
B.B. Жолобов*, Е.И. Тарновскии
МОДЕЛИРОВАНИЕ НЕУСТАНОВИВШИХСЯ ТЕЧЕНИЙ УГЛЕВОДОРОДНЫХ СМЕСЕЙ
В ТРУБОПРОВОДАХ
"Управление магистральных нефтепроводов Центральной Сибири, Томский государственный педагогический университет "Томский политехнический университет
Трубопроводный транспорт, являясь доминирующим видом перемещения энергоносителей, играет особую роль для нашей страны. Современная наука и техника не нашли альтернативы трубопроводам, поэтому их значение в ближайшие десятилетия будет лишь возрастать. Эксплуатация трубопроводных систем порождает большое число вопросов, правильно ответить на которые можно, опираясь только на научные проработки. Для теоретического исследования течения в трубах, как правило, используются модели механики сплошных сред. Физические законы сохранения, примененные к элементарному объему транспортируемого флюида, приводят к определяющей системе дифференциальных уравнений, гибкость которой обеспечивается введением ряда параметров, идентифицируемых впоследствии по диспетчерской информации. Неустановившиеся напорные и безнапорные течения в гидравлическом приближении описываются системами дифференциальных уравнений в частных производных гиперболического типа [1].
Уравнения в частных производных этого типа достаточно хорошо изучены и разработан ряд эффективных численных методов для их решения. Привлекательной стороной такого подхода является универсальность и возможность адаптации к реальным условиям. Традиционно напорные и без-
напорные течения в трубопроводном транспорте рассматривались изолированно. Это, по-видимому, связано с тем, что сопряженное решение задачи о возникновении и схлопывании паровых или парогазовых полостей слабо поддается аналитическому анализу. Имеются лишь отдельные разрозненные публикации, в которых изучается поведение изолированных газопаровых полостей в ограниченных объемах жидкости. В случае существенного изменения давления при движении углеводородных смесей и при неустановившихся процессах возникают фазовые переходы, приводящие к образованию паровых и парогазовых объемов различного масштаба (от изолированных пузырьков в объеме жидкости до протяженных участков с неполным заполнением поперечного сечения трубы). Полости, возникающие при смене режимов перекачки, при восстановлении исходного режима могут исчезать (схлопываться) в том случае, если они заполнены только паровыми компонентами углеводородной жидкости. Различие в механизмах выделения и растворения газовых компонентов приводит к тому, что для растворения газовых компонентов необходимо уже настолько больше времени, что реализуются режимы течения с неполным заполнением трубы при давлениях, существенно превышающих давление насыщенных паров жидкости.
Для прогнозирования процессов, происходящих при траспортировке нефти в этом случае уже необходимо знать неравновесный компонентный состав фаз. Первичный компонентный состав газопаровой фазы, образовавшейся в волне разрежения, зависит от компонентного состава нефти и может быть определен путем рассмотрения фазового равновесия нефти как многокомпонентной углеводородной смеси. Последующее трансформирование компонентного состава газопаровой фазы необходимо рассматривать в рамках неравновесной схемы фазовых переходов, учитывающей предысторию процесса. Таким образом, одним из важных этапов построения математической модели, описывающей движение многокомпонентной смеси, является конкретизация функций, задающих интенсивность межфазного переноса компонентов.
Применение принципов термодинамики необратимых процессов [2] позволяет получать некоторые феноменологические линейные соотношения для кинетики фазовых переходов [3]. При этом феноменологические коэффициенты, входящие в эти соотношения, зависят от свойств фаз, структуры смеси, геометрии области движения и должны определяться из эксперимента или детального анализа микродвижений в окрестности поверхностей фазового перехода. Молекулярно-кинетичес-кий анализ процессов уноса и осаждения молекул пара на единицу площади межфазной границы приводит к формуле Герца-Кнудсена-Ленгмюра [4, 5, 6] для предельной скорости фазовых превращений. Использование этого уравнения в инженерных целях требует знания коэффициента аккомодации для составляющих углеводородную смесь компонентов. Относительно этого коэффициента не только не существует теории, пригодной для предсказания его величины, но и отсутствуют простые способы ее экспериментального определения. Имеющиеся в настоящее время данные о коэффициенте аккомодации для различных веществ очень противоречивы. В связи с этим имеет смысл рассмотрение упрощенной квазиравновесной схемы на межфазной границе, предполагающей, что, несмотря на неравновесность в объеме фаз, на самой границе выполняются условия термодинамического равновесия.
Дальнейшим упрощением является принятие более сильного допущения в виде квазиравновесной схемы в объеме непрерывных фаз, предполагающей выполнение условий термодинамического равновесия в объеме смеси одной из фаз с каплями или пузырьками другой. Это упрощение предполагает, что сначала каждая из фаз автономно подстраивается к изменению внешних условий, быстро трансформируя свой компонентный состав в направлении локального термодинамического равновесия внутри занимаемого ею объема. На втором (существенно более длительном) этапе
устанавливается истинное термодинамическое равновесие между сосуществующими фазами. В реальных условиях транспортировки истинное фазовое равновесие недостижимо. Кинетика межфазного массообмена в элементарном объеме определяется микроскопическими движениями компонентов смеси и структурой межфазной границы. Можно ожидать, что эти движения должны описываться соответствующим решением уравнений многокомпонентной диффузии. Даже в простейших случаях это можно осуществить лишь при условии правильного задания граничных условий.
Учитывая сложную форму межфазных границ смеси в реальной углеводородной смеси, которую практически невозможно описать в математическом виде, не удается решить и задачу определения интенсивности межфазного массообмена каким-либо точным на микроуровне способом. В связи с этим полагается, что доминирующими факторами, определяющими интенсивность массообмена г —» /, являются скорости образования промежуточной фазы г—»у"в дисперсном виде и скорости перехода её в непрерывную форму г" —» у (осаждение капель, слияние пузырьков). По сравнению с характерным временем процесса транспортировки, временем образования и расслоения промежуточных фаз можно пренебречь и принять макроскопически изотермическую и локально равновесную по фазовому переходу г -» у* модель течения смеси газа, паров и капель нефти. Интенсивность массообмена между непрерывной и дискретной формами определяет в этом случае эффективную скорость массообмена в модели раздельного движения фаз. Принимая гравитационный механизм осаждения (всплытия) капель (пузырьков), введем характерные времена в{ процесса объединения дисперсной и непрерывной частей фаз. Предполагая, что в момент времени / в объеме непрерывных частей фаз существуют только те дисперсные включения, которые образовались во временном промежутке {-в., получим следующие соотношения для интенсивностей межфазного массообмена:
(1)
где Уу - вертикальная скорость движения дисперсных включений от у'-й фазы к г'-й фазе, которая определяется следующими обыкновенными дифференциальными уравнениями
V»
й1
-М/
Г по
„О
1
Ч р-,1
о.
У
Х] Й „о
С1и Рг
7 = 0,1),
где % ~ длина линии контакта непрерывных фаз, а о - объемная концентрация пузырьков пара в жидкости, а] - объемная концентрация капель
в газопаровой фазе, р 9 - истинная плотность соответствующей фазы в непрерывной форме.
В данной работе предполагается, что в трубопроводе имеются напорные участки, заполненные жидкостью целиком, а также участки, заполненные лишь частично, давление на которых равно давлению в парогазовой фазе. В отличие от работы [7] рассматривается случай, когда давление в парогазовой фазе может отличаться от давления насыщенных паров жидкосги. Ограничимся рассмотрением изотермических течений вязкой сжимаемой жидкости. Систему уравнений для средних гидродинамических величин, выражающую балансовый закон сохранения массы и количества движения фаз, можно представить следующим образом [1]:
5 С . 3 17
ОТ ОХ
р = (1 — а )р°, (/ Ф ] = ОД),
1 У 1
дж
рр{
дх
+
^РМ^ + -г р$хРхи! = ~р1
ОТ дх
дГ.
+ (ж -ж ) —'- -рр, дх
- \4idXi X,
(2)
я^р + р.^.,
здесь р- - средняя по сечению кажущаяся плотность г'-й фазы (включающей дисперсную часть у'-й фазы - газопаровая фаза с каплями жидкости или жидкая фаза с пузырьками пара); ^ - часть поперечного сечения трубы, занятая соответствующей фазой (г'-й фазой в непрерывной форме, включая дисперсную частьу-й фазы); и{ - средняя по поперечному сечению трубы скорость г'-й фазы; (5г - коэффициент, характеризующий поправку на количество движения потока жидкости, которая протекает в единицу времени через поперечное сечение с площадью вследствие неравномерности распределения скорости и плотности по сечению (/?,=1 в случае идеальной жидкости, при турбулентном движении /3 = 1.03-ь1.1 при неустановившемся движении /Забудет переменной величиной, зависящей от характера распределения скоростей в сечениях трубы); гг— расстояние по вертикали вверх от некоторой горизонтальной плоскости до данной точки г'-й фазы; 71 х - среднее значение величины 7Г,-по периметру сечения соответствующего потока; г0г- - величина силы трения, действующая на единицу площади боковой поверхности трубы; р -давление, которое предполагается одинаковым
на границе раздела обеих фаз; %г- - периметр сечения потока г'-й фазы.
Предполагая, что давление по сечению жидкой фазы распределяется по гидростатическому
закону: р,= р0 - Р^2 ~ <),
где г2- расстояние по вертикали вверх от некоторой горизонтальной плоскости до свободной поверхности жидкости; g - ускорение свободного падения; р] - давление в парогазовой фазе, которое можно предполагать слабо меняющимся по поперечному сечению. Тогда
я: =^ + р2§г2 =рх
ж
ср„
Ж
ср =кХх =Р 1
(3)
так как является функцией только х. Если ввести в рассмотрение глубину потока жидкости Л, а через гД обозначить вертикальную координату дна, то
г2 =2д+И
с!х
и
дг*2 _ с12Д д дх с1х дх
Д
\2
(Зх
где г'г - гидравлический уклон на безнапорном участке течения (уклон дна). Часть периметра сечения потока, который примыкает к твердым стенкам (смоченный периметр), обозначим через ¿¡. Тогда гидравлический радиус потока опре-
делится выражением гГ
I.
. Объемная концен-
трация газовой фазы и смоченный периметр связаны с глубиной потока жидкости в трубопроводе следующими соотношениями:
V-
>0
)Р = (Л - Я)^ - (И - я/ +
п2 • Н-Я _2 ж
+ Л агсвт-+ Я —
Я 2
Ь1 =2Лагсвш
Ь-Я К
+ лЛ,
Я»
К +.Р, =Р = кИ\ ¿0+Ц = 2кК,
(4)
га-Робинсона [10], имеющее для индивидуального вещества следующий вид:
Р
Е*Т
(о - Ь) п(о + Ь) + Ъ(ь - Ь)
(7)
здесь - внутренний радиус трубы; объемное содержание парогазовой фазы в непрерывной форме, включающей капли; площадь поперечного сечения трубопровода.
Второе уравнение системы (2) с учетом (3) можно представить в виде
3 д
РхЩЪ + -Г" Ъ =
сж дх ох
а = аса{Т), а
0.4572Я*2ГС2
, 0.073Л Г „
1 + т
Г Т\°-5
т = 0.37464 + 1.54226ПУ - 0.26992т'
сЬс Эх
<к
\2
Д
сЬс
(5)
- \гыс1хг -иур^и).
Сила трения, действующая на единицу площади боковой поверхности трубы, может быть представлена следующим образом:
1 2 1 Т0. = р.и. Д.
рг
Ке^е,-,
I2 '
(6)
р{щ4гГ
где Ле,- = ■
Л-
2г.
г.
I
Здесь £г-- относительная шероховатость, е -средний размер выступов шероховатостей.
Для замыкания системы (1-6) воспользуемся уравнениями состояния фаз. Кроме этого, положим, что объемные концентрации дисперсных включений равны их мольным долям при локальном термодинамическом равновесии фаз. Определение составов и количественного соотношения равновесных паровой и жидкой фаз (при заданном давлении, температуре и общем составе смеси) осуществляется на основе классических положений термодинамики фазового равновесия многокомпонентных смесей, неоднократно рассмотренных в соответствующей литературе.
В настоящее время в качестве уравнений состояния углеводородных смесей широко используются уравнения состояния ван-дер-ваальсовс-кого типа. Наиболее удачные из большого количества имеющихся модификаций обсуждаются в [9]. В данной работе используется уравнение Пен-
0.424 Г, 10.2/0 Ъ ш =----0.994, =
(1 -тк) к Тс
где Т- абсолютная температура, К; о - мольный объем, м3/моль; Я* - универсальная газовая постоянная; со- ацентрический фактор; Ть- температура кипения; индекс (с) приписывается параметрам в критической точке.
Для расчета смесей применяются комбинационные правила:
где % - условное обозначение мольной доли £-го компонента парогазовой или жидкой фазы; с^- коэффициенты парного взаимодействия, значения которых для ряда индивидуальных углеводородов приведены в работе [11].
Выражение для летучести А:-го компонента, необходимое для записи условий равновесия, имеет вид [9]
Ч=\п(г}о)-)ф-В) +
ЪЛг-1)
-А
2Е
П,{\-сы){акап)
0.5
а
1п-
г — 2.4Ш] 2-0.418 I
А--
а р
т
¥4^
, в =
Ъ р
ГУ!
Я*Т
, (к ф я =1,2,
здесь суммирование ведется по числу компонентов N.
Коэффициент сжимаемости г определяется из кубического уравнения
г2 - (1 - В)г2 +{А- ЪВ2 - 2В)г -
-(АВ-В2 -В3) = 0.
При рассмотрении парогазовой фазы смеси выбирается максимальный действительный корень, а при рассмотрении жидкой фазы - минимальный действительный корень. Из условий термодинамического равновесия, материального баланса между фазами и нормировки следует замкнутая система уравнений для расчета локального равновесного состава фаз:
4
о
4' ^М У
И«
ax)xkl,
Чkl = Ук2 V-а2 У+ а2^2' 2 XU У Id =
летучести к-го компонента в рав-
г о л здесь У И и ь
новесной парогазожидкостной смеси с капельной и пузырьковой структурой соответственно; хЫ'УИ'Г11а - мольные доли к-го компонента г'-й фазы в газопаровой, жидкой составляющей фаз и смеси непрерывных и дисперсных составляющих соответственно; а1 ~ мольные доли общего количества вещества г'-й фазы в непрерывной форме; £ - суммирование по индексу к. Текущая величина Г)к1 вычисляется с помощью массовых концентраций С, определяемых из уравнений сохранения массы компонент
д д
от ах
Ух
„, „Отг хк2Мк2
Щр\ v\
My
М-,
(9)
здесь Мк1 - молекулярный вес компонентов г'-й фазы, М. - молекулярный вес вещества дисперсных включений соответствующей фазы. При записи зависимостей (9) принято, что объемное содержание фаз в дисперсной форме равно их текущей молярной концентрации в соответствующей несущей фазе.
Коэффициент динамической вязкосги фаз газовых и жидких углеводородных смесей, входящее в соотношение (6), при давлениях больших 1МПа с учетом введенных ранее обозначений рассчитывается с помощью соотношений
п =Р01 +1-08 * 1(Г4 [ехр(1 Л39рпр1 )-— ехр(-1.111/>^58 )] <§.,
где - коэффициент динамической вязкости при атмосферном давлении и температуре
/% = 3.4 *1(Г4; рпр1,Тпр1 - приведенные плотность и температура
npi
Pi
'пщп
1 I /пкр1
' "Р1 т .тш
L npi
mpi
1
пкрг iV
м 1J
N 7=1
т • = УпТ ••
>мпкр1 ppij'
Prucpi ^ ^nKpiTnKpi/VnKpi,
N
znKpi = HljZij.Mi
j=i
J У'
N i=i
Зависимость площади сечения трубы от давления для упругого материала стенок, подчиняющегося закону Гука, выражается следующим соотношением:
F = Fn
1 +
Dy (Р-Ро)
Е
где Бу - внутренний диаметр трубы, § - толщина стенки трубы, Е - модуль упругости материала стенок трубы, ^ - площадь поперечного сечения трубы при отсутствии избыточного внутреннего давления.
Далее будем полагать, что до момента времени ( = 0 движение отсутствовало или было стационарным. По отношению к массовому расходу С> и давлению р имеем следующие начальные условия:
г = 0, в=в-во=0, р=р-р0=0.
Граничные условия зависят от вида возмущений на границах течения. Во многих задачах трубопроводного транспорта на одном конце трубопровода известно давление как функция времени. К другому концу трубопровода присоединен насосный агрегат, изменяющий расход жидкости по известному закону. При этом насосный агрегат может быть присоединен к трубопроводу непосредственно или отделен от него устройством для регулирования расхода или уменьшения колебаний давления. Граничные условия в этих случаях могут быть записаны в виде [1]
зё
х = 0,р= (p(t), x = l,Q+g
дх
где (p(t) и y/(t) - известные функции, равные нулю при /=0; g - постоянная величина. Другой тип граничных условий задается в виде соотношения, связывающего расход и давление:
д = т(?)^Ахр + А2,
где т(/)~ заданная функция времени; Ах, Аг - некоторые постоянные. Это условие имеет место, например, в сечении, где нарушается герметичность трубопровода, при открытии (закрытии) задвижки и в некоторых других случаях.
Для численного расчета течений нефти или нефтепродуктов как на полностью заполненных участках трубопровода, так и на заполненных лишь частично используется аналог метода распада произвольного разрыва С.К. Годунова [12]. Рассматривая метод распада разрыва как формальную процедуру, приводящую к устойчивой разностной схеме, разностные уравнения для системы (2, 5, 9) представляются в виде, аналогичном случаю идеального газа. «Большие» значения определяются из условий распада произвольного разрыва в предположении отсутствия силы трения. Распадные значения для каждой из фаз в отдельности определяются из распада разрыва в идеальном газе, аппроксимирующем соотношения (7). При этом используется не истинная плотность фаз, а кажущаяся плотность, учитывающая наличие дисперсных включений. Уравнения для определения параметров на границах ячеек в этом случае совпадают с уравнениями для расчета распада разрыва в равновесном запыленном газе [13]. Особенность заключается в определении величин ^ на границе расчетной ячейки. Для жидкой фазы оказывается эффективным использование направленной схемы [13] типа
п+0.5 (У2 )л+0.5 +(^2>Л+1.5 - О/
л+1.5. (^г)п+0.5 +(^2)и+1.5 -О-
Для газопаровой фазы считается, что распадные значения устанавливаются в минимальном проходном сечении )п+1 = тш.{(^ )п+0.5'(^1 )п+1.5 }•
Выбор шага по времени осуществляется таким образом, чтобы ни в каком месте разностной сетки не успевало произойти взаимодействие волн, образующихся в результате распада разрыва на границах соседних ячеек.
Таким образом, задача о движении многокомпонентной углеводородной жидкости при нестационарных режимах трубопроводной транспортировки в математическом плане сводится к решению системы уравнений (2, 5, 9) с соответствующими начальными и граничными условиями. При этом, в отличие от [7], нет необходимости вводить различные системы уравнений для описания движения на напорных и безнапорных участках.
Численный расчет движения и распределения многокомпонентной углеводородной смеси, состоящей из модельного состава, представленного в таблице, в трубопроводе при внутренном давлении р = 0.6 МПа, ¿=293 Кис граничными условиями соответствующие свободному истечению с одного торца трубопровода с р = 0.3 МПа и поступлению смеси при р = 0.6 МПа осуществлялся на основе принципа установления равновесного состояния по всей длине трубопровода.
Для всех расчетов профиль трубопровода 5=1 м выбирался с геодезическими координатами, представленными схематически (рис. 1).
Компоненты Пк
СН4 0.1
с?нЙ 0.0482
С}нн 0.0293
С¡Н12 0.3684
С7н,6 0.1323
С10Н22 0.2252
Н28 0.0966
Ь*5 м
Рис. 1. Профиль трубопровода
Шм
Рис. 2. Распределение фазового давления в трубопроводе при истечении многокомпонентной смеси
Рис. 4. Распределение скорости течения жидкой фазы в трубопроводе
S,%
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Ь*5 м
Рис. 5. Положение границы раздела жидкой и паровой фаз в различные моменты времени
На рис. 2 приведены результаты расчета распре- и жидкой фаз соответственно по длине трубо-
деления давления в фазах вдоль профиля трубы провода. На рис. 5 - положение границы разде-
на выделенных участках. На рис. 3, 4 представ- ла жидкой и паровой фаз в различные моменты
лено распределение скорости течения паровой времени.
Литература
1. Гинзбург И.П. Прикладная гидрогазодинамика. Л., 1958.
2. Базаров И.П. Термодинамика. М.,1983.
3. Нигматулин Р.И. Динамика многофазных сред. Ч. 1. М., 1987.
4. Hertz Н. Ober die Verdunstung der Flüssigkeiten, insbesondere des Quecksilbers im luftleeren Räume // Der Physik und Chemie. 1882. Bd. 17,177.
5. Knudsen M. Die maximale verdampfungsgeschvindigkeit des Quecksilbers // Annalen der Phusik. 1915. Bd. 47. № 13.
6. Langmuir I. The Dissociation of Hydrogen into atoms // The Journal of the American Chemical Society. 1915. V. 37. № 3.
7. Агавин A.A., Карасевич A.M., Сухарев М.Г. и др. Трубопроводные системы энергетики: модели, приложения, информационные технологии. М., 2000.
8. Рид Р. и др. Свойства газов и жидкостей. Л., 1982.
9. Баталии О.Ю. и др. Фазовые равновесия в системах природных углеводородов. М., 1992.
10. Peng D.Y., Robinson D.B. A new two-constant equation of Stele // Ind. End. Chem. Fundam. 1976. V. 15. № 1.
11. Qerlich L., Plocker U., Prausnitz J.M., Knapp H. Equation-of-state methods for computing Phase equilibria and endhalpis // International Chemical Enginering. 1981. № 1.
12. Годунов C.K. Разностный метод численного расчета разрывных решений гидродинамики // Математический сборник. 1959. Т. 47. № 3.
13. Матвеев С.К. Некоторые аспекты применения метода Годунова к решению нестационарных задач газовой динамики // Уч. зап. Ленингр. ун-та. Сер.: Газодинамика и теплообмен. 1977. № 5.
14. Многомерная и многокомпонентная фильтрация / С.Я. Закиров, Б.Е. Сомов, В.Я. Гордон. М., 1988.
