Численный анализ процесса парообразования при кипении аварийного пролива сжиженного газа Текст научной статьи по специальности «Физика»

Е. В. Старовойтова

аспирант Казанского государственного технологического университета, г. Казань, Татарстан

А. Д. Галеев

канд. техн. наук, доцент Казанского государственного технологического университета, г. Казань, Татарстан

С. И. Поникаров

д-р техн. наук, профессор Казанского государственного технологического университета, г. Казань, Татарстан

УДК 536.24:614.83

ЧИСЛЕННЫЙ АНАЛИЗ ПРОЦЕССА ПАРООБРАЗОВАНИЯ ПРИ КИПЕНИИ АВАРИЙНОГО ПРОЛИВА СЖИЖЕННОГО ГАЗА

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

Ключевые слова: аварийный пролив; сжиженный газ; кипение; численное моделирование.

Введение

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

Можно выделить следующие процессы, определяющие динамику развития аварии с СГ:

• истечение СГ из технологических систем в случае их повреждения или внезапного полного разрушения;

• мгновенное вскипание СГ;

• разлив жидкой фазы на подстилающую поверхность и кипение СГ;

• рассеивание паровых облаков в граничном слое атмосферы;

• сгорание облака топливно-воздушной смеси. Недостаточное знание этих процессов не только

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

Одним из ключевых моментов при оценке последствий возможных аварий на объектах хранения, транспортирования и переработки сжиженных

© Старовойтова Е. В., Галеев А. Д., Поникаров С. И., 2011

газов является определение количества вещества, поступающего в окружающее пространство с поверхности пролива при кипении СГ, с учетом следующих факторов:

• свойств подстилающей поверхности;

• преобладающих атмосферных условий (температура воздуха, скорость ветра и т. д.);

• изменения температуры жидкой фазы.

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

Математическая модель парообразования сжиженного газа

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

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

24

0869-7493 ПОЖАРОВЗРЫВОБЕЗОПАСНаСТЬ 2011 ТОМ 20 №2

В данной статье предлагается модель, описывающая процесс парообразования СГ при аварийном проливе. При разработке модели были сделаны следующие допущения:

• теплоприток от подстилающей поверхности определяется теплопередачей по твердой фазе;

• не учитывается фазовый переход влаги в подстилающем слое.

Интенсивность парообразования W (кг/(м2-с)) рассчитывалась следующим образом:

W =

аИ

ц

W

пРи qs + qa ^ ; пРи qs + qa < ит>

(1)

где qa — тепловой поток из атмосферы, Вт/м2; qs — тепловой поток от подстилающей поверхности к жидкости, Вт/м2; АН — теплота парообразования, Дж/кг;

—диффузионный поток, определяемый с помощью пристеночных функций из условия, что мольная доля пара на межфазной границе У* « 1 (в модели принималось У* = 0,98), кг/(м2-с);

— массовый поток пара при диффузионном испарении, кг/(м2-с).

Первое условие в формуле (1) соответствует режиму кипения, при котором температура жидкости остается постоянной; второе условие — режиму диффузионного испарения, при котором изменение температуры жидкости рассчитывается по формуле

¿Т1 = qs + qa ~ »У АИ dх С Р1 ш1

(2)

где Т1 — температура жидкой фазы, К; х — время, с;

СР1 — удельная теплоемкость жидкости, Дж/(кг-К);

т1 — масса жидкости, отнесенная к единице площади поверхности пролива, кг/м2. Мольная доля пара на межфазной границе в режиме диффузионного испарения рассчитывалась по формуле

У* =

- в*

(Т)

Рп

(3)

где У* — мольная доля бутана на границе раздела фаз;

Рш (Т) — давление насыщенных паров бутана при температуре Т, Па; Р0 — давление в окружающей среде, Па. Массовый поток пара с поверхности пролива определялся на основе стандартных функций стенки [1] с учетом поправки на стефановский поток [2]:

= к&

(СУ ~ СР ) ри*

с + ;

(4)

С + = <1 Бс У + при у + < уС;

_Бс{(и++ РС) при у+> уС;

+ Ри* УР

у =-;

= - 1п(Еу + ) - АВ; к

(5)

(6) (7)

АВ =

0 при г + < 2,25;

-

+

г01 - 2,25

л

+ Сгг 0+

87,75

V У

X бш{0,4258 (1п г + - 0,811)} при 2,25 < г + < 90;

-1п(1 + Сгг0+) при г0+ > 90; -

(8)

0

р г 01 и*

РС = 9,24

Бс,

3/4

-1

К5 =

[1 + 0,28 е-°'0018с/8с'

1п [1(1 - У*)]

(9) (10) (11)

где К8 — коэффициент, учитывающий влияние сте-фановского потока на интенсивность испарения; С* — концентрация пара на границе раздела фаз, кг/кг;

СР — концентрация пара в пристеночном узле расчетной сетки, кг/кг;

р — плотность паровоздушной смеси, кг/м3; и* — скорость трения; и* = (х*/р)0'5;

напряжение трения на стенке, Па;

Бс, Бс — молекулярное и турбулентное числа Шмидта соответственно; Бс = ц/(рВт); Бт—коэффициент молекулярной диффузии, м2/с; уР — расстояние по нормали от поверхности испарения (стенки) до соседнего узла расчетной сетки, м;

ц — коэффициент молекулярной динамической вязкости, кг/(м-с);

уС+ — безразмерное расстояние от стенки, определяемое в точке пересечения линейного и логарифмического закона изменения концентрации у стенки;

к — константа Кармана, равная 0,41; Е — константа в логарифмическом законе стенки для скорости, равная 9,1; Сг — коэффициент;

г01 — высота шероховатости поверхности пролива, м;

г+ — безразмерная высота шероховатости поверхности пролива.

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

0869-7493 ПОЖАРОВЗРЫВОБЕЗОПАСНаСТЬ 2011 ТОМ 20 №2

25

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

В программе FLUENT для учета влияния шероховатости подстилающей поверхности следует задать два параметра — высоту шероховатости z0 и константу шероховатости Cz. Соотношение, связывающее аэродинамическую шероховатость y0 c z0, имеет вид [3]:

, _ 9,793уо

(12)

Значение аэродинамической шероховатости для поверхности пролива у0 = 0,001 м [4], для поверхности, прилегающей к проливу (равнинная местность, сельскохозяйственные угодья), у0 = 0,025 м [4]. Таким образом, при расчете принимали С2 =1, высоту шероховатости для поверхности пролива хы = = 0,0098 м и для прилегающей территории — г0 = = 0,245 м.

При использовании пристеночных функций требуется, чтобы расстояние от стенки до соседнего узла расчетной сетки удовлетворяло условию уР > г0 (рис. 1). Учитывая данное обстоятельство, вертикальный размер пристенных ячеек со стороны атмосферы задавали равным 0,5 м.

Тепловой поток, подводимый от подстилающей поверхности к жидкой фазе (д5 = /5у)у = 0),

определялся из численного решения трехмерного нестационарного уравнения теплопроводности для твердого слоя:

C P 9Ts К

Csp S -ft _ К

^ 2TS

dx 2

g 2ts

dy 2

^

dz 2

(13)

где Т8(х, у, £)—распределение температур в подстилающем слое;

С8, р,, ^ — теплоемкость, плотность и коэффициент теплопроводности твердого материала соответственно.

При решении уравнения (12) на твердой поверхности, прилегающей к проливу, задавалось условие

?;=т.

Тепловой поток из атмосферы да вычислялся с помощью пристеночных функций [1].

Для учета дополнительного нагрева воздуха вследствие конденсации водяного пара использовалась функция источникового члена: • в уравнении энергии:

Qv _

V(T)СУ раНу ; Ат '

(14)

А # р

Ур W

'////vs/////-'//

Рис. 1. Расположение пристеночного узла расчетной сетки

• в уравнении переноса компонента (паров воды) ¥(Т)СУ рАЯу _

5 _ —

Ат

(15)

где у(Т) — функция, определяющая долю сконденсированной влаги;

'0, Т > Тл;

273,15 < Т < Тё; (16)

1, Т < 273,15 К;

V (T) _

CV CVd Су

Т — температура паровоздушной среды, К; Тй — температура точки росы при данном парциальном давлении паров воды в паровоздушной смеси, К;

- концентрация паров воды в воздухе, кг/кг; уЛ — концентрация водяного пара, при кото-

Су C

рой он достигает насыщения при температуре Т,

кг/кг;

Ат — шаг по времени, с;

АНГ — теплота парообразования воды, Дж/кг.

Результаты расчетов

Исходные данные для расчетов принимались следующие: сжиженный газ — бутан; начальная температура жидкости 272,65 К; АН = 385600 Дж/кг; СР1 = 2300 ДжДкг-К); Бп = 9,44-10-6 м2/с; Бс, = 0,7; размеры аварийного пролива 20 х 20 м; температура окружающей среды 311 К; скорость ветра варьировалась и имела значения 1, 2, 3 м/с на высоте 10 м; начальная толщина пролива 0,05 м; шероховатость поверхности проливахы = 0,01 м; шероховатость поверхности территории, прилегающей к проливу, = 0,25 м. Свойства твердого материала: = 1,28 Вт/(м К); С3 =1130 Дж/(кгК); р, = 2300 кг/м3; АНГ = = 2,45-106 Дж/кг.

Проводилось сравнение результатов расчета удельной массы испарившегося бутана тисп, полученных с использованием разработанной модели и формулы [5]:

ти

(Ts - Tl)

аН

2KJ— na

5,h/Re Кa

(17)

ISSN 0869-7493 ООЖАРОВЗРЫВОБЕЗООАСООСТЬ 2011 ТОМ 20 №2

о

, кг/м2

Ти К

12 10 8 6 4 2

0 500

,сп> кг/м2

1000 1500 2000

х, с

14 12 10 8 6 4 2

О 500 1000 1500 2000 х, с

'"исп. кг/м2

18 16 14 12 10 8 6 4 2

500

1000

1500 2000

-с, с

Рис. 2. Зависимость удельной массы тисп испарившегося бутана от времени х при скорости ветра: а — 1 м/с; б — 2 м/с; в — 3 м/с; 1 — по численному методу; 2 — по формуле (17)

где а — эффективный коэффициент температуропроводности материала, на поверхность которого произошел разлив вещества, м2/с; Яе — число Рейнольдса; Яе = Ш/\а; и — скорость воздушного потока, м/с; d — характерный размер пролива СУГ, м; \а — кинематическая вязкость воздуха, м2/с; ~ка — коэффициент теплопроводности воздуха, Вт/(м-К).

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

Значения тисп, рассчитанные по численной модели, при скорости ветра 1 м/с ниже, а при скоростях 2 и 3 м/с выше соответствующих значений, определенных по формуле (17). Это объясняется тем, что в формуле (17) не учитывается переход к диффу-

1000 1500 2000

т,с

Рис. 3. Изменение температуры сжиженного бутана Т во времени х: 1,2,3 — скорость ветра соответственно 1; 2 и 3 м/с

Чи, Вт/м2 12000 10000 8000 6000 4000 2000

3 2 1

0

500 1000 1500 2000 т, с

?а. Вт/м2

Рис. 4. Зависимость тепловых потоков qs и qa от време -ни х: а — от подстилающей поверхности; б — из атмосферы; 1,2,3 — скорость ветра соответственно 1;2и3м/с

зионному режиму испарения, который имеет место на поздних стадиях процесса парообразования при высоких скоростях ветрового потока, когда диффузионный поток превышает тепловые потоки от подстилающей поверхности и воздуха. При этом различие в значениях за рассматриваемый промежуток времени составляет: при скорости ветра 1м/с — не более 10 %; 2 м/с — не более 13 %; 3 м/с — не более 23 %. При прочих равных условиях удельная масса испарившегося бутана тем выше, чем больше скорость ветра. Так, при скорости ветра 3 м/с величина тисп за рассматриваемый промежуток времени примерно на 38 % выше, чем при скорости ветра 1 м/с.

На рис. 3 представлены зависимости температуры жидкости от времени, полученные по численной

0869-7493 ПОЖАРОВЗРЫВОБЕЗОПАСНОСТЬ 2011 ТОМ 20 №2

27

модели. При низкой скорости ветра (1 м/с) температура в течение длительного времени не изменяется, так как суммарный теплоприток из окружающей среды превосходит теплоотвод вследствие диффузионного испарения. При скоростях ветра 2 и 3 м/с температура начинает снижаться через 100 и 60 с соответственно.

На рис. 4 показано изменение тепловых потоков от подстилающей поверхности и атмосферы к жидкой фазе во времени т. В режиме кипения при рассматриваемых скоростях ветра значения д8 совпадают. Снижение температуры жидкости при переходе к диффузионному режиму испарения обуславливает более высокие значения д8 при скоростях ветра 2 и 3 м/с, чем при 1 м/с. Увеличение дасо временем связано с повышением температурного градиента со стороны атмосферы.

Доля теплопритока из атмосферы в суммарном теплопритоке за рассматриваемый промежуток времени составляет не более 4,5 % при скорости ветра 1 м/с, 11,6 % — 2 м/с, 17,2 % — 3 м/с.

Заключение

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

СПИСОК ЛИТЕРАТУРЫ

1. Fluent Inc. Fluent 6.1. User's Guide, Lebanon, 2003.

2. Померанцев В. В., Арефьев К. М., Ахмедов Д. Б. Основы практической теории горения. — Л.: Энергоатомиздат, 1986. — 312 с.

3. Blocken В., Stathopoulos Т., Carmeliet J. CFD simulation of the atmospheric boundary layer: wall function problems // Atmospheric Environment. — 2007. — Vol. 41. — Issue 2. — P. 238-252.

4. РД 03-26-2007. Методические указания по оценке последствий аварийных выбросов опасных веществ : утв. приказом Ростехнадзора 14 декабря 2007 г. № 859 // Экологические ведомости. — 2008. — № 4.

5. ГОСТР 12.3.047-98. ССБТ. Пожарная безопасность технологических процессов. Общие требования. Методы контроля. — Введ. 2000-01-01. — М. : Изд-во стандартов, 1998.

Материал поступил в редакцию 15 ноября 2010 г.

Электронный адрес авторов: starovojtova@inbox.ru.

2!

Ш ISSN 0869-7493 ООЖАРОВЗРЫВОБЕЗООАСНОСТЬ 2011 ТОМ 20 №2