Компьютерное моделирование процессов распространения в атмосфере газовых выбросов при авариях на газопроводах

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

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

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

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

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

Модель основана на решении краевой задачи для одномерного уравнения турбулентной диффузии следующего вида:

где C(t,x) - одномерная концентрация загрязнения в точке х и в момент времени t, Q - масса аварийного выброса, который первоначально сосредоточен в начале координат, К - коэффициент турбулентной диффузии, ?(*) - 8- функция Дирака.

Задача (3.2.2) имеет известное аналитическое решение в виде функции гауссовского типа:

В трехмерном случае уравнение турбулентной диффузии в приближении с постоянными коэффициентами турбулентности также имеет аналитическое решение. Краевая задача в этом случае имеет вид:

где Кх, Ку, К_ - коэффициенты турбулентной диффузии вдоль координатных осей.

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

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

В частном случае заметного преобладания ветрового переноса над переносом за счет турбулентной диффузии уравнение (3.2.6) принимает упрощенный вид:

Рассмотрим уравнение одномерного переноса аварийного облака (3.2.6) и выполним его компьютерное моделированием средствами MatLab. Для компьютерного моделирования (3.2.6) будем использовать метод конечных разностей. При этом область моделирования разбивается на конечное число узлов с шагом Ах, в которых ищется решение уравнения. Компьютерное моделирование временной динамики также осуществляется с постоянным временным шагом At на временных слоях. При моделировании адвективного члена будем использовать явную вычислительную схему, которая получила название схемы разностей против потока (РПП):

Решение уравнения (3.2.8) получается в виде итерационной формулы:

где Ки = — - параметр Куранта - Фридрихса - Леви, D-К- безраз-

Ах Ах

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

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

где Sk - величина ошибки на к-п итерации.

Для устойчивости вычислительной схемы необходимо, чтобы модуль множителя перехода был меньше единицы:

Решение уравнения (3.2.6) можно разложить в ряд Фурье, поэтому устойчивость вычислительной схемы (3.2.9) исследуем с помощью произвольной Фурье-моды следующего вида:

где a(t) - амплитуда Фурье-моды, р=— - волновое число, X - длина вол-

Л

ны Фурье-моды, i - мнимая единица.

Подставим (3.2.12) в уравнение (3.2.9) и сократим на е'рх. В результате получим:

Из уравнения (3.2.13) следует выражение для множителя перехода:

где введено обозначение р = рАх.

Преобразуя экспоненты в (3.2.14) по формулам Эйлера, получим:

Условие устойчивости (3.2.11) принимает следующий вид:

Выражение, стоящее в левой части неравенства (3.2.16), принимает максимальные значения при cos/? = ±1. В обеих случаях sin/? = 0.

При cos/?= 1 из (3.2.16) получим следующую величину для множителя перехода |g|2 = 1. При этом вычислительная схема попадает на границу устойчивости.

При cos/? = -l из (3.2.16) получаем следующее условие устойчивости:

Из (3.2.17) следует:

Сократив (3.2.18) на положительный множитель K„+2D, получим:

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

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

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

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

Далее приведены тексты программ обеих модулей.

%TURB_SNR1 распространение загрязнений от мгновенного источника

N=10; % размер сетки по оси ОХ

М - 12; % кол-во итерационных шагов

dx = 1; % шаг сетки по оси ОХ

dt = 0.4; % временной шаг модели

и = 2; % скорость ветра

К = 0.2; % коэффициент турбулентной диффузии

Ки - u*dt/dx; % Ки = 0.8 - параметр (Куранта-Фридрихса-Леви)

D - K*dt/dx/dx; % D - 0.08 - безразмерный коэффициент турбулентной диффузии

% численная схема устойчивая: Ки + 2D - 0.96 < 1 С = turbdiffl (N, М, dx, dt,Ku,D);

% построение графиков п = 1:N;

V = С(2,п); VI = С(3,п); V2 = C(4,n); V3 = С(5,п); V4 = С(6,п); figure('NumberTitle', 'off, 'Name', 'Расчет турбулентной диффузии'); hPlot=plot(n, V,n, VI,п, V2,n, V3,n, V4); title('Pac4em турбулентной диффузии методом РПП); grid on;

xlabel(Расстояние вдоль оси ОХ);

у1аЬе1('Концентрация );

text(4.7,1.25, 'leftarrow2-H итерация);

text(5.6,l, 'leftarrow3-H итерация);

text(6.4,0.85, 'leftarrow4-я итерация);

return

%TURBDIFF1 численное решение одномерного уравнения переноса при постоянной

% скорости ветра и постоянном коэффициенте турбулентной диффузии

% Деформация фронта загрязнения.

% Граничные условия 1-го рода, function С = turbdiffl (N,M,dx,dt,Ku,D)

N - размер сетки по оси ОХ М- кол-во итерационных шагов %dx- шаг сетки по оси ОХ %dt- временной шаг модели

%Ки - параметр устойчивости (Куранта-Фридрихса-Леви)

%oD- безразмерный коэффициент турбулентной диффузии С = zeros(M,N); % инициализация массива

С(1,:) = 0; % инициализация начального распределения загрязнения

С(1,4) = 2; % начальная концентрация в аварийном облаке

С(1,5) = 2;

С(:,1) - 0; % граничные условия

C(:,N) = 0; j = 2:N-1;

% главный итерационный блок модели for к= 1:М-1

C(k+l,j) = C(k,j) - Ku*(C(k,j) - C(k,j-1)) + ...

D*(C(kJ+l) + C(kJ-l) - 2*C(kj)); end return

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

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

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

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

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