Численное моделирование загрязнения акватории Р. Днепр при аварийной утечке аммиака из аммиакопровода «Тольятти Одесса» Текст научной статьи по специальности «Физика»

УДК 519.6 : 504.453

А. Н. ПШИНЬКО, Л. В. АМЕЛИНА, Н. Н. БЕЛЯЕВ (ДИИТ)

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ЗАГРЯЗНЕНИЯ АКВАТОРИИ р. ДНЕПР ПРИ АВАРИЙНОЙ УТЕЧКЕ АММИАКА ИЗ АММИАКОПРОВОДА «ТОЛЬЯТТИ - ОДЕССА»

На 6a3i розроблено! 2D чисельно! моделi виконано розрахунок забруднення aKBaTopii' р. Дшпро при ава-рiйному витоку aмiaку. Модель базуеться на чисельному штегруванш рiвняння конвективно-дифузiйного переносу домiшки та моделi потенцiйного руху. Наводяться результати обчислювального експерименту.

Ключовi слова: аваршний витiк aмiaку, забруднення акватори, ашакопровод, конвективно-дифузiйний перенос домiшки, двовимiрнa чисельна модель

На базе построенной 2D численной модели выполнен расчет загрязнения р. Днепр при аварийной утечке аммиака. Модель основывается на численном интегрировании уравнения конвективно-диффузионного переноса примеси и модели потенциального течения. Приводятся результаты вычислительного эксперимента.

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

The 2D numerical model was developed and used to simulate river pollution after accident on the ammonia pipe over Dnipro River. The model is based on the numerical integration of the K-gradient transport model and potential flow. The results of numerical experiment are presented.

Keywords: ammonia emergency leakage, water pollution, ammonia pipe, convection-diffusion transfer of impurity, two-dimensional numerical model

Введение

В настоящее время одной из актуальных экологических проблем является прогноз последствий аварий на крупных химически опасных объектах. Одним из таких объектов на территории Днепропетровской области есть аммиакопровод (АП) «Тольятти - Одесса». Эта транспортная система открыто пересекает р. Днепр ниже г. Днепропетровска (рис. 1). Очевидно, что в случае аварии на этом участке АП произойдет масштабное загрязнение окружающей среды.

Рис. 1. Вид аммиакопровода «Тольятти - Одесса» при пересечении р. Днепр

Одной из важных задач при оценке последствий такой аварии является прогноз уровня загрязнения акватории р. Днепр. Данный про-

гноз может быть выполнен только на основе метода математического моделирования. Учитывая, что в настоящее время предъявляются высокие требования к качеству прогнозной информации, что проявляется в первую очередь в необходимости наиболее полного учета в прогнозной модели факторов, влияющих на процесс переноса загрязнителя, применение для рассматриваемого прогноза эмпирических или аналитических моделей - является нецелесообразным [1, 8]. Цель данной работы - разработка численной модели загрязнения акватории реки при попадании в нее сжиженного аммиака из АП. Основной интерес при моделировании на базе построенной модели представляет оценка глубины проникновения струи аммиака в воду. Для прогноза в работе используется 2Б модель. Как известно двухмерные модели продолжают оставаться эффективным инструментом исследования многопараметрического процесса массопереноса [10].

Математическая модель процесса

Будем считать, что в случае аварии на АП образовалось отверстие в трубопроводе и струя сжиженного аммиака начала поступать в акваторию (рис. 2).

Для расчета распространения загрязнителя в акватории реки используется подход, примененный ранее в ГИАПе [6], когда тепловые

© Пшинько А. Н., Амелина Л. В., Беляев Н. Н., 2010

процессы при моделировании загрязнения водной среды в случае попадания в нее аммиака не рассматривались. В данной работе для описания процесса транспорта загрязнителя в акватории (профильная задача) используется осред-ненное по ширине загрязненной струи уравнение переноса, выражающее в дифференциальном виде закон сохранения массы для загрязнителя [4, 2, 8]:

dC duC dvC

dt dx

dy

= div (| gradC),

(1)

где С - концентрация загрязнителя в воде; и, V, - компоненты вектора скорости потока; ц = = (цх , цу) - коэффициент турбулентной диффузии; ^ - время; х, у - декартовы координаты. Ось «у» направлена вертикально вверх, ось «х» - вдоль потока.

C| г

= C

E '

На выходной границе расчетной области, в дискретной модели ставится «циклическое» (мягкое) граничное условие вида

C (i +1, j) = C (i, j),

где i, j - номер соответствующей разностной ячейки.

В построенной численной модели учитывается неравномерный профиль скорости течения воды в реке. Расчет профиля скорости водного потока осуществлялся по следующей зависимости [1, с. 80-81]:

= u0 • >/1 - Pz2

где u0 - скорость течения у поверхности;

(H - y)

z = -

H

- относительная глубина;

Рис. 2. Расчетная схема аварийной утечки аммиака из аммиакопровода «Тольятти - Одесса»

Постановка краевых условий для данного уравнения рассмотрена в работах [2, 4]. Отметим, что дно реки и свободная поверхность являются граничными линиями тока. В построенной численной модели на этих границах реализуется граничное условие вида

дС=о,

дп

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

На входной границе полагается, что концентрация загрязнителя равна нулю. На том участке верхней границы, где в акваторию поступает загрязнитель - струя аммиака (рис. 2) ставится граничное условие вида:

Н - глубина потока по вертикали; у - текущая координата (начало координат у дна);

Р - параметр.

Величина Р рассчитывалась по формуле

Р = 0,0222С - 0,000197С2,

где С - коэффициент Шези.

Данная зависимость рекомендуется при 60 < С < 90 .

Для расчета скорости течения у поверхности используется зависимость [8]

С

V к 111___V

С -1 '

где Vср - средняя скорость течения.

Коэффициенты диффузии аппроксимировались следующими зависимостями

ц х = 1,5 • В V;

I z =

g • H •Vx

48 • C

где CE - известное значение концентрации загрязнителя.

где В - ширина реки.

Из поврежденного АП струя аммиака выходит из отверстия с определенной скоростью, изменяющейся со временем в связи с падением давления в поврежденном АП. В первые секунды истечения эта скорость будет максимальной, и ее значение будем использовать для прогноза глубины проникновения струи аммиака в акваторию. Скорость истечения аммиака из поврежденной трубы (полагаем, что с такой скоростью или меньшей струя аммиака взаимодействует со свободной поверхностью) рассчитывается по зависимости [7]:

Vв

= ф-л/2g н

где ф - коэффициент скорости;

Н - напор.

Особая сложность состоит в задании величины коэффициента скорости, т.к. форма отверстия (и его размеры) в поврежденном АП заранее неизвестны. Однако очевидно, что это будет отверстие неправильной формы, возможно - щель, а значит, для такого случая коэффициент скорости будет минимальной величиной.

Отметим, что в работе [7] скорость истечения поврежденного АП составляла величину порядка 42 м/с (при среднем давлении в аммиа-копроводе 70 ати), а коэффициент скорости был равен 0,59. Учитывая, что за переходом АП через р. Днепр находится станция подкачки, можно считать, что давление в АП на открытом участке через р. Днепр ниже этой средней величины. Поэтому в данной работе расчет выполнен для минимального давления в АП -15 ати.

Для применения модели (1) необходимо рассчитать поле скорости потока, формирующееся при взаимодействии входящей струи аммиака и течения в р. Днепр. Для расчета поля скорости при таком взаимодействии применяется модель невязкой жидкости, и моделирующим уравнением является уравнение для потенциала скорости [3]:

д2 Р д2 Р

дх2 ду

= 0.

(2)

= др = дР_ дх ду

Метод решения

Для численного интегрирования уравнения переноса используется попеременно-треугольная разностная схема [2]. Разностные соотношения в операторном виде записываются так [2]:

, 1

- на первом шаге расщепления к = п +— :

1 - j + 1(L+Ck + L+yCk) =

At 2

= \(MXxCk + MXxCk + M+yC" + MyCn);

, 1

- на втором шаге расщепления k = n + ;

c=n+—: 4

Ck-CC 1( At 2 ~

- + ^ (L-Ck + LLyCk) =

Для данного уравнения ставятся следующие граничные условия [2, 3]:

дР

- на дне и свободной поверхности — = 0 , где

дп

п - единичный вектор внешней нормали к границе;

- на входной границе (и на участке верхней границы, где имеет место втекание аммиака

дР

в акваторию) — = Vn , где Vn - известное

дп

значение скорости потока (скорости втекания);

- на выходной границе расчетной области Р = Р * (х = const, у) + const (условие Дирихле).

В рамках рассмотренной модели компоненты вектора скорости рассчитываются по зависимостям [3]:

= !(M„Ck + M+Cc + M-yCk + M+yCC);

- на третьем шаге расщепления k = п + ^ ;

c=п+—: 2

Ck cc i

11 - 1 +1 (L+Ck + L-Ck) = At 2 x y

=4( m„c с + Myck + M-yCk + M+yCc)

- на четвертом шаге расщепления k = п +1;

3

c=п+—:

4

Ck - Cc 1

-Л-L + i( LZCk + L+Ck) =

At 2 x y

1 —i

4

= ^ (MXxCk + M+Cc + MyCc + M++yCk).

В данных выражениях использованы следующие обозначения разностных операторов:

1

du+C « u+1,-u<+Cnw = L+cn+l-

Щ+1 .

dy

дг Ax

дu C - f^n+1 n+1 ^ ui+1,jCi+1,j -uijCij

дx Ax

ду+C v + Cn+1 - v+Cn+1 „ 4 j+1^ij ji, j-1

дy Ay

3v~C v- Cn+1 - v"C^1 „ 4j+1^i,j+1 vij4j

= L-C ,

= L+C

n+1.

Ay

= L-C

n+1 .

осуществляется по явной формуле - методу бегущего счета [9].

После определения величины потенциала скорости компоненты вектора скорости на гранях ячеек рассчитываются по соотношениям:

uj =

vj =

P - P Ax

P - P

ri,j ri, j-1

Ay

д , dC _ Ci+1, j - Cij —(Д x—) ~ Д x— 2—-

M x ^ ' JX .7

дх dx Ax

-ДXC] ~C-1,j =M-C +M+xC;

Ax

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

dC

C - C

—(д —) «д "iJ+1 -д/^ д/ Ay2

-Д = M-yC +M++C-

Ax2

yy

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

Для численного интегрирования уравнения (2) используется метод установления решения по времени. Поэтому численно интегрируется уравнение эволюционного вида

P

It

д2 P д2 P

дт2 дТ

(3)

где ^ - фиктивное время.

При / ^ да решение уравнения (3) будет стремится к «установлению», т.е. к решению уравнения (2).

Для численного интегрирования уравнения (3) используется неявная схема суммарной аппроксимации [9]. В этом случае разностные уравнения на каждом дробном шаге имеют вид:

n+— P 2 - Pn i, j_ij

At

1

n+—

- P .2

i, j

+ P 2

i n+—

Pn+1 - P 2

i, j i, j

At

Ax2

pn+1 pn+1 P+1, j - Pi, j

1 1 n+— n+— _P 2 + P 2

4, j j-1

Ay

2

Ax2

pn+1 pn+1 Pi, j+1 - Pi, j

Ay2

Данная схема является неявной, но расчет значения потенциала скорости р ^ в каждой разностной ячейке на каждом шаге расщепления

Практическая реализация модели

На основе построенной численной модели создан код на алгоритмическом языке FORTRAN. Разработанный код построен на модульном принципе, что позволяет легко настроить его на решение широкого круга задач данного класса. Данный код был использован для расчета загрязнения акватории р. Днепр и, в первую очередь, для прогноза глубины проникновения струи аммиака в реку. При проведении вычислительного эксперимента использовались данные, предоставленные проектно-изыскательским институтом «Днепрогипровод-хоз» (г. Днепропетровск), а также данные литературных источников - разработок ГИАПа [7]. Расчет выполнен при таких параметрах: средняя скорость течения в р. Днепр на рассматриваемом участке 0,07 м/с; коэффициент скорости 0,5; скорость входа струи аммиака в воду -10 м/с, давление в АП - 15 ати; глубина реки 29 м; ширина реки - 900 м; концентрация загрязнителя в струе 611 кг/м3; температура аммиака в АП 20 °С; длина расчетной области 85 м. Т.к. АП располагается на высоте 16 м от уровня воды, то за счет торможения воздухом произойдет уменьшение скорости входящей струи. Поэтому расчет также выполнен для меньшей скорости входящей в воду р. Днепр струи аммиака - 7 м/с.

Результаты вычислительного эксперимента показаны на рис. 3 - 7. На данных рисунках представлены изолинии концентрации аммиака в реке для различных моментов времени (здесь представлено распределение концентрации в процентах от концентрации во входящем потоке). Представленные рисунки позволяют выявить динамику формирования зоны загрязне-

ния в акватории в первые секунды после аварии на АП. Отчетливо видно, что формирующаяся в акватории реки зона загрязнения разворачивается за счет взаимодействия двух потоков по аналогии сноса ветром факела загрязнителя, выходящего из устья трубы. На представленных рисунках хорошо видно, что в силу неравномерного течения в реке зона загрязнения существенно втягивается в сторону течения вблизи поверхности. Уменьшение скорости входящей струи аммиака на 30 % (с 10 м/с до 7 м/с) не привело к существенному уменьшению размеров зоны загрязнения. Так, например, для момента времени t = 10 с глубина зоны загрязнения для скорости входа струи 10 м/с составляет величину порядка 22 м, а при скорости входа 7 м/с - порядка 18 м. Формирование данной зоны связано с достаточно большой кинетической энергией входящего в акваторию реки потока аммиака. Таким образом, при аварии произойдет загрязнение акватории р. Днепр практически на большую часть ее глубины (подчеркнем, что расчет сделан для минимального давления в аммиакопроводе на момент аварии).

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

1.112Е*01 о 1.104Е+01 г 1.972Е+00 а 1.900Е+00 1 1.В20Е«00 и 1.75БЕ+00 а 1.Б84Е+00 t 1.612Е+00 в 1.540Е+00 1.4БВЕ+00 1.39БЕ+00 у 1.324Е+00 ].252ЕЛ ). 180Е+00 ). 10БЕ+00 1.3Б0Е-01

ЯШ ~ 25

^10

40

новения загрязнителя в реку. Расчеты показали, что глубина проникновения загрязнителя составит порядка 22 м. Учитывая, что средняя глубина реки на рассматриваемом участке акватории составляет порядка 29...30 м, очевидно, что в случае аварийного повреждения ам-миакопровода произойдет интенсивное загрязнение акватории.

Рис. 3. Зона загрязнения акватории реки, t = 2 с (скорость струи аммиака 10 м/с)

Выводы

На основе построенной численной модели проведен вычислительный эксперимент по оценке уровня загрязнения акватории р. Днепр в случае аварийной утечки аммиака при повреждении открытого участка аммиакопровода «Тольятти - Одесса», расположенного над рекой. Основной интерес при моделировании представлял вопрос прогноза глубины проник-

Рис. 4. Зона загрязнения акватории реки, t = 5 с (скорость струи аммиака 10 м/с)

Рис. 5. Зона загрязнения акватории реки, t = 10 с (скорость струи аммиака 10 м/с)

Рис. 6. Зона загрязнения акватории реки, t = 10 с (скорость струи аммиака 7 м/с)

Дальнейшее развитие рассмотренного в работе направления необходимо проводить в направлении создания плановой модели рассеивания загрязнителя в русле совместно с моделью рассеивания загрязнителя в атмосфере для комплексной прогнозной оценки последствий аварии на аммиакопроводе «Тольятти -Одесса».

0.690Е+00 ____

0.170Е+01 г.ппг(1 ¡пл+я х 0.033Е*02

Рис. 7. Зона загрязнения акватории реки, t = 15 с (скорость струи аммиака 10 м/с)

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Караушев, А. В. Сборник задач по гидравлике [Текст]. - Ч. 2 / А. В. Караушев, Н. А. Панчу-рин. - Л.: Речной транспорт, 1957. - 196 с.

2. Численное моделирование распространения загрязнения в окружающей среде [Текст] / М. З. Згуровский [и др.] - К.: Наук. думка, 1997. - 368 с.

3. Лойцянский, Л. Г. Механика жидкости и газа [Текст] / Л. Г. Лойцянский. - М.: Наука, 1978. -735 с.

4. Марчук, Г. И. Математическое моделирование в проблеме окружающей среды [Текст] / Г. И. Марчук. - М.: Наука, 1982. - 320 с.

5. Методические основы оценки и регламентирования антропогенного влияния на качество поверхностных вод [Текст] / под ред. А. В. Ка-раушева. - Л.: Гидрометеоиздат, 1987. - 285 с.

6. Проникновение аммиака в грунт и его распространение в грунтовых водах [Текст] / С. В. Не-рпин [и др.] // Магистральные трубопроводы: Труды ГИАП. - Вып. 51. - М., 1978. - С.43-66.

7. Несвижский, Ф. А. Расчет длины зоны поражения [Текст] / Ф. А. Несвижский, А. П. Потанин, Н. С. Эвенчик // Магистральные аммиакопро-воды: Труды ГИАП. - Вып. 51. - М., 1978. -С. 34-43.

8. Основы прогнозирования качества поверхностных вод [Текст]. - М.: Наука, 1982. - 181 с.

9. Самарский, А. А. Теория разностных схем [Текст] / А. А. Самарский. - М.: Наука, 1983. -616 с.

10. Tedeschi, G. Study of vertical transport of marine aerosol using an unsteady 2D model [Text] / G. Tedeschi // Conf. Abstracts of 31st NATO/SPS Int'l Techn. Meeting on Air Pollution Modelling and Its Application (27 Sept. - 01 Oct. 2010, Torino, Italy). - 2010. - P. 4.9.

Поступила в редколлегию 18.06.2010.

Принята к печати 25.06.2010.