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

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

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

где x, у, z - пространственные координаты в декартовой прямоугольной системе координат, u,v,w - компоненты скорости ветрового потока, р - отклонение давления от гидростатического, р - средняя плотность воздуха для области решения, К - коэффициент турбулентной вязкости.

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

Для решения поставленной задачи можно использовать алгоритм MAC. Рассмотрим трехмерный вариант этого алгоритма. В алгоритме MAC используется итерационная процедура, состоящая из двух этапов. На первом этапе осуществляется решение первых трех уравнений системы (3.5.1), а на втором этапе решается уравнение для давления, которое можно представить в следующем виде:

_ di dv dw

где D — —+— + — - дивергенция скорости.

дх, dy dz

Уравнение (3.5.2) представляет собой уравнение Пуассона, для решения которого разработано множество алгоритмов. Большинство этих алгоритмов представляют модификации четырех основных методов: частично неявный итерационный метод последовательной верхней релаксации, частично неявный итерационный метод чередующихся направлений, метод циклической редукции и метод рядов Фурье. Как отмечалось выше, успех применения того или иного метода обусловлен в основном правильной постановкой граничных условий. Рассмотрим метод последовательной верхней релаксации.

Граничные условия для уравнения давления (3.5.2) ставятся в виде нулевых условий Неймана на верхней и выходной границах. На входной границе задается граничное значение давления. На нижней границе и стенках элемента рельефа задается градиент давления, который определяется из уравнений (3.5.1), записанных для точек границы. Так для вертикальной наветренной границы граничное условие для давления имеет следующий вид:

где Рф - отклонение давления от гидростатического в граничном узле, i,

j, к - индексы узлов расчетной сетки, Ах, Ay, Az - размеры ячейки расчетной сетки по х, у и z соответственно, g - индекс граничного узла.

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

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

где С, = -Spp - правая часть уравнения Пуассона, Д = Дх/Ду, Д, = Axl Az, со - параметр релаксации.

Рассмотрим численные результаты компьютерного моделирования распространения аварийного облака с учетом обтекания элементов рельефа. В качестве модельного элемента рельефа рассмотрим орографический барьер в форме параллелепипеда. Область моделирования размером 400x400x200 м разделим на элементарные ячейки с помощью конечно-разностной сетки размером: 40x40x20 ячеек. Таким образом, расчетная сетка будет представлять собой набор кубических ячеек со стороной Юм. Введем орографический барьер размером 100x70x60 м в форме параллелепипеда с координатами: (21 - 30), (21 - 27), (1 - 6). Алгоритм включает две итерационные процедуры. Основная итерационная процедура, которая решает первые три уравнения системы (3.5.1), и вспомогательная итерационная процедура, реализующая метод последовательной верхней релаксации для расчета полей давления. Вспомогательная итерационная процедура осуществляется на каждом из шагов основной итерационной процедуры. При этом в алгоритме для расчета давления используется параметр верхней релаксации со = 1,2.

Критерием сходимости алгоритма верхней последовательной релаксации является максимальное рассогласование:

где «1 - номер итерации, el - оценка точности сходимости итерационной процедуры.

Для основной итерационной процедуры используется аналогичный критерий, который для и - компоненты скорости потока имеет следующий вид:

Временной шаг для основной итерационной процедуры равен: At = 2 с. На рис.3.5.1 показан график Е(п) сходимости основной итерационной процедуры. В расчетах использованы следующие значения параметров: К - 1 м2/с, р = 1 кг/мЗ, начальное и граничное условия для скорости (u,v,w) - (2,0,0) м/с, начальное и граничное условия для давления р = 1 Па (такая величина давления выбрана вследствие того, что система уравнений (2) не зависит от абсолютных величин давления, а лишь от градиента давления).

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

График сходимости основной итерационной процедуры при компьютерном моделировании распространения аварийного облака с учетом обтекания препятствия, п - номер итерации

Рис.3.5.1. График сходимости основной итерационной процедуры при компьютерном моделировании распространения аварийного облака с учетом обтекания препятствия, п - номер итерации

На рис.3.5.2 показаны траектории облака, как для потока, полученного по уравнениям (3.5.1), так и для потенциального приближения. Для расчета потенциального течения, аналогично Хино , примем, что потенциал течения через элемент рельефа может быть представлен в виде потенциала двойного слоя по следующей формуле:

где: h(x,y) - высота поверхности, D - область течения, U - скорость набе-

I 2 2 2

тающего потока, r = y](x-a) + (у~Р) +z расстояние от точки (а,/3,0) до точки (x,y,z).

В приближении точечного диполя из (3.5.7) могут быть получены следующие упрощенные формулы для расчета потенциального течения:

где hk и Sk - высота и площадь к-ого препятствия, х ,у - координаты

U/C У) К

центра масс к-ото препятствия, гк = +(т~Тм)' +(Z_Z0J" •

Траектория аварийного облака при обтекании элемента рельефа в форме параллелепипеда и при отсутствии мелкомасштабных пульсаций

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

Для идентификации модели в качестве параметров настройки можно использовать следующие два параметра: lk =hkSk и zok=z'-z, где z' - смещенная

вертикальная координата.

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

где ри - автокорреляционная функция Лагранжа для i-компоненты скорости ветра, S. - случайная компонента скорости ветра, - стандартное отклонение для /-компоненты скорости ветра, Тц - Лагранжев временной масштаб автокорреляции, и' - мелкомасштабные пульсации скорости ветра, At - интервал автокорреляции.

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

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

Траектории аварийных облаков при обтекании препятствия с учетом мелкомасштабных пульсаций

Рис.3.5.3. Траектории аварийных облаков при обтекании препятствия с учетом мелкомасштабных пульсаций: а - потенциальное приближение, б - результат применения гидродинамической модели (3.5.1)

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

(рис.3.5.2) поглощается пульсационной компонентой (рис.3.5.3), что в конечном итоге должно приводить к удовлетворительному согласию между результатами расчетов полей атмосферного переноса газовых облаков, выполненных в рамках этих приближений. Изолинии полей распространения газового облака, рассчитанные по 20 реализациям траекторий, как для случая потенциального течений, так и для случая течений Рейнольдса, показаны на рис.3.5.4.

Изолинии концентрации аварийного облака при обтекании препятствия (интенсивность выброса - 0,8 кг/с)

Рис.3.5.4. Изолинии концентрации аварийного облака при обтекании препятствия (интенсивность выброса - 0,8 кг/с): а - потенциальное приближение, б - гидродинамическая модель; 1 -1,0-10"5; 2 - 2,2-10’5; 3 - 4,4-10'5; 4 -1,0-10'4; 5 -1,8-10’4; 6 - 2,2-Ю'4 в кг-м'3. Пунктирная линия изображает начальный подъем струи на эффективную высоту

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

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

где N - число узлов расчетной сетки, q^k, q™ - объемная концентрация

газовых выбросов в узле ijk соответственно для потенциального течения и течения Рейнольдса.

Расчет критерия s осуществлялся для различных высот источника газовых выбросов. Результаты представлены в виде графика на рис.3.5.5.

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

Отличие двухмерных полей аварийного облака для потенциального течения и гидродинамической модели

Рис.3.5.5. Отличие двухмерных полей аварийного облака для потенциального течения и гидродинамической модели: h - эффективная высота подъема газовой струи

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

 
Посмотреть оригинал
< Пред   СОДЕРЖАНИЕ   ОРИГИНАЛ     След >