Прогнозирование последствий аварийных залповых выбросов сжиженных газов (стр. 2 )

.

(10)

Изменение массы сжиженного газа в жидкой фазе рассчитывалось из уравнения:

.

(11)

Уравнения (10) и (11) были интегрированы в пакет FLUENT с помощью пользовательских функций (UDF).

Массовый поток газа с поверхности пролива определялся из уравнений:

,	(12)
	(13)
,	(14)
,	(15)
,	(16)

где Kstef – коэффициент, учитывающий влияние на интенсивность испарения стефановского потока; Yg,s– концентрация пара на границе раздела фаз, кг/кг; Yg,P– концентрация пара в пристеночном узле расчетной сетки, кг/кг; Сμ – константа; kр – турбулентная кинетическая энергия в пристеночном узле, м2/с2; Y+ – безразмерная концентрация; Sc, Sct – молекулярное и турбулентное числа Шмидта соответственно; yС+ – безразмерное расстояние от стенки, определяемое в точке пересечения линейного и логарифмического закона изменения концентрации у стенки; yР – расстояние по нормали от поверхности испарения (стенки) до соседнего узла расчетной сетки, м; m – коэффициент молекулярной динамической вязкости, кг/(м∙с); k – константа Кармана, равная 0,41; E – константа в логарифмическом законе стенки для скорости, равная 9,1.

Поправка на стефановский поток (16) была введена в стандартные функции стенки, встроенные в пакет FLUENT, с помощью UDF.

Тепловой поток, подводимый от грунта к жидкой фазе, qgrd = lgrd(¶Tgrd/¶y)y=0 определялся из численного решения трехмерного нестационарного уравнения теплопроводности для твердого подстилающего слоя:

,

(17)

где CР,grd, rgrd, lgrd, Tgrd – соответственно теплоемкость, плотность, коэффициент теплопроводности, температура подстилающей поверхности.

Граничные условия для уравнения (17) задавались следующие (см. рис.2):

· на твердой поверхности, прилегающей к проливу Tgrd(t,x,0,z)=Tliq,

· на нижней и боковых границах подстилающего слоя ¶Tgrd/¶n=0, где n – нормаль к поверхности, ограничивающей расчетную область.

Распределение температур в начальный момент времени в подстилающем слое: Tgrd (0, x, y, z)=Ta, где Ta – температура окружающего воздуха, К.

Тепловой поток из атмосферы qа вычислялся с помощью пристеночных функций. Тепловые потоки qs, qp и qar определялись с помощью простых формул, приведенных в литературе.

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

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

Пользовательские функции использовались для учета дополнительного нагрева воздуха вследствие конденсации водяного пара в холодном аэрозольном облаке путем включения источниковых членов Qv в уравнении энергии и Sv в уравнении переноса компонента (паров воды):

; ,

(18)

где DHν – теплота парообразования воды, Дж/кг; ψ – разность между концентрацией паров воды в воздухе и концентрацией насыщения, кг/кг; ρ – плотность паровоздушной среды, кг/м3.

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

Разработанная модель парообразования тестировалась как по результатам проведенного экспериментального исследования процесса кипения жидкого азота на непроницаемой поверхности, описанного в Главе 2, так и по экспериментальным значениям, опубликованными в открытой печати (рис. 3). В целом наблюдается хорошее согласование экспериментальных данных с расчетными данными, вычисленными по представленной модели. Для проверки адекватности модели парообразования на стадии испарения с поверхности пролива однокомпонентной жидкости использовались экспериментальные данные по испарению сжиженного бутана, опубликованные в открытой печати. Результаты сравнения представлены в табл.1.

Для сравнения в табл. 1 представлены значения, вычисленные по формуле:

,

(19)

где М – молекулярный вес вещества, кг/моль; U – скорость движения воздуха, м/с; Pg,s – давление насыщенных паров при начальной температуре, мм. рт. ст.

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

Рисунок 3 – Сравнение результатов расчетов по разработанной модели с экспериментальными данными: - полученными автором (сжиженный азот);

представленными в работах: - со сжиженным азотом*; - со сжиженным метаном**

Таблица 1

Экспериментальные и расчетные значения интенсивности испарения

№ эксп	Эксперимен-тальное значение интенсивности испарения, г/(м2·c)***	Интенсивность испарения, рассчитанная по модели, г/(м2·c)	Отклонение расчетных значений от эксперимен-тальных, в %	Интенсив-ность испарения, рассчитан-ная по ур.(19), г/(м2·c)	Отклонение расчетных значений по ур.(19) от эксперимен-тальных, в %
1	3,7	3,2	-15,6	0,72	-413,9
2	2,7	3,5	+22,9	0,81	-233
3	2,3	2,0	-15	0,63	-265
4	3,9	2,7	-44,4	1,3	-200
5	5,5	4,5	-22,2	1,1	-400
6	5,6	6,4	+12,5	1,1	-409
7	7,2	6,6	-9,1	1,2	-500

Для подтверждения адекватности модели распространения токсичного газа в атмосфере проведен сравнительный анализ расчетов с данными аварии с аммиаком, происшедшей 13 июля 1973 г. в Потчефструме (ЮАР). На рис. 4 представлена схема расположения погибших и пострадавших на месте аварии.

Границы зоны токсического поражения аммиаком оценивались по значению показателя токсодозы Dg(x,y,z):

,

(20)

где tэксп — время экспозиции, с; N — количество шагов по времени; Dt — шаг по времени, с; Yg — концентрация опасного вещества в воздухе, кг/кг.

Функция вычисления токсодозы была реализована на языке программирования С (Си) и интегрирована в пакет FLUENT.

Как видно из рис. 4, все погибшие находятся в предсказанной зоне смертельного поражения. Общая зона поражения складывается за счет прохождения двух облаков: первичного и вторичного, последнее образуется при испарении пролива. От первичного облака формируется основной круг поражения. Вторичному облаку соответствует выступ в направлении ветра на контуре зоны поражения. Форма и размеры зон летального поражения при 10 и 30 минутах отличаются незначительно. Отсюда следует, что смертельная интоксикация аммиаком обусловлена, в большей степени, распространением первичного облака в течение первых минут после аварии.

Рисунок 4 – Схема расположения пострадавших, в том числе погибших непосредственно на месте аварии в Потчефструме (ЮАР, 13.07.73).

Обозначения: Внутренний контур – граница зоны смертельного поражения при времени экспозиции tэксп= 10 мин, внешний контур – граница зоны смертельного поражения при tэксп= 30 мин: 1 – разрушенный резервуар; 2 – местоположение людей, найденных мертвыми; 3 – местоположение людей, впоследствии умерших; 4 – местоположение людей во время аварии

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

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

Таблица 2

Удельная масса испарившегося аммиака из пролива (t=1800с)

Параметр	Значения
j=0%	j=50%	j=100%
1 м/с	2 м/с	3 м/с	1 м/с	2 м/с	3 м/с	1 м/с	2 м/с	3 м/с
mисп, кг/м2	6,6	7,7	8,6	6,9	7,8	8,7	7,4	8,0	8,8

Примечание. j − относительная влажность воздуха.